%\documentclass[11pt,a4paper, oneside]{article}
\documentclass[nojss]{jss}
%\documentclass[article]{jss}
\usepackage{thumbpdf,lmodern}
\usepackage{Sweave}
\usepackage{amsmath,amsfonts,enumerate}
\usepackage{bm}
\usepackage{mathrsfs}
\usepackage{bbm}
\usepackage{tabularx}
\usepackage{multirow}
\usepackage{xcolor}
\definecolor{Green}{RGB}{0,150,0}


\newcolumntype{C}{>{\centering\arraybackslash}X}
%%\VignetteIndexEntry{mvord: An R Package for Fitting Multivariate Ordinal Regression Models}
%%\VignetteDepends{mvord}
%% need no \usepackage{Sweave.sty}
\author{Rainer Hirk\\
WU Wirtschaftsuniversit\"at\\
Wien
\And Kurt Hornik\\
WU Wirtschaftsuniversit\"at\\
Wien
\And Laura Vana\\
WU Wirtschaftsuniversit\"at\\
Wien}
\Plainauthor{Rainer Hirk, Kurt Hornik and Laura Vana}

\Address{
Rainer Hirk\\
Department of Finance, Accounting and Statistics\\
Institute for Statistics and Mathematics\\
WU Wirtschaftsuniversit\"at Wien\\
1020 Vienna, Austria\\
E-mail: \email{rainer.hirk@wu.ac.at}\\ \ \\
Kurt Hornik\\
Department of Finance, Accounting and Statistics\\
Institute for Statistics and Mathematics\\
WU Wirtschaftsuniversit\"at Wien\\
1020 Vienna, Austria\\
E-mail: \email{kurt.hornik@wu.ac.at}
\\ \ \\
Laura Vana\\
Department of Finance, Accounting and Statistics\\
Institute for Statistics and Mathematics\\
WU Wirtschaftsuniversit\"at Wien\\
1020 Vienna, Austria\\
E-mail: \email{laura.vana@wu.ac.at}
}

\title{\pkg{mvord}: An \proglang{R} Package for Fitting Multivariate Ordinal Regression Models}
\Plaintitle{mvord: An R Package for Fitting Multivariate Ordinal Regression Models}

\Abstract{The \proglang{R} package \pkg{mvord} implements composite
  likelihood estimation in the class of multivariate ordinal
  regression models with probit and a logit link.
  % error structures
  A flexible modeling framework for multiple ordinal measurements on
  the same subject is set up, which takes into consideration the
  dependence among the multiple observations by employing different
  error structures.
  % covariate dependent errors
  Heterogeneity in the error structure across the subjects can be
  accounted for by the package, which allows for covariate dependent
  error structures.
  % constraints
  In addition, regression coefficients and threshold parameters are
  varying across the multiple response dimensions in the default
  implementation. However, constraints can be defined by the user if a
  reduction of the parameter space is desired.
}

\Keywords{composite likelihood, cross-sectional data, longitudinal
  data, multivariate ordinal regression model, \proglang{R}}
\Plainkeywords{composite likelihood, cross-sectional data,
  longitudinal data, multivariate ordinal regression model, R}
<<echo = FALSE, results = hide>>=
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
@
\begin{document}
\SweaveOpts{concordance=TRUE}
<<echo=FALSE>>=
library("mvord")
data("data_cr_panel")
data("data_cr_mvord")
data("data_cr_mvord2")
cache <- TRUE
@
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Introduction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\sloppy
% Intro
% ------------
The analysis of ordinal data is an important task in various areas of
research. One of the most common settings is the modeling of
preferences or opinions (on a scale from, say, \emph{poor} to
\emph{very good} or \emph{strongly disagree} to \emph{strongly
  agree}).  The scenarios involved range from psychology (e.g.,
aptitude and personality testing), marketing (e.g., consumer
preferences research) and economics and finance (e.g., credit risk
assessment for sovereigns or firms) to information retrieval (where
documents are ranked by the user according to their relevance) and
medical sciences (e.g., modelling of pain severity or cancer stages).

Most of these applications deal with correlated ordinal data, as
typically multiple ordinal measurements or outcomes are
available for a collection of subjects or objects (e.g., interviewees
answering different questions, different raters assigning credit
ratings to a firm, pain levels being recorded for patients repeatedly
over a period of time, etc.).
% Approach
% ------------
In such a multivariate setting, models which can deal with the
correlation in the ordinal outcomes are desired. One possibility is to
employ a \emph{multivariate ordinal regression model} where the
marginal distribution of the subject errors is assumed to be
multivariate. Other options are the inclusion of random-effects in the
ordinal regression model and conditional models
\citep{fahrmeir2001multivariate}.

% R packages
% ------------
Several ordinal regression models can be employed for the analysis of
ordinal data, with cumulative link models being the most popular
ones \cite[see][]{Tutz12, ordinal}. Other approaches include
continuation-ratio or adjacent-category models \cite[see][]{Agresti02, Agresti10}.  Different packages to analyze and
model ordinal data are available in \proglang{R} \citep{RCoreTeam}:
For univariate ordinal regression models with fixed effects the
function \code{polr()} of the \pkg{MASS} package \citep{MASSpkg}, the
function \code{clm()} of the \pkg{ordinal} package \citep{ordinalR}, which supports scale effects as well as nominal effects,
and the function \code{vglm()} of the \pkg{VGAM} package
\citep{VGAMpkg} are available. Another package which accounts for
heteroskedasticity is package \pkg{oglmx}
\citep{oglmx}. The package \pkg{ordinalNet} \citep{ordinalNetpkg}
offers tools for model selection by using an elastic net penalty. Moreover, ordinal logistic models can be fitted by the functions \code{lms()} and  \code{orm()} in package \pkg{rms} \citep{rmspkg} and the \code{MCMCoprobit()} function in package \pkg{MCMCpack} \citep{MCMCpackpkg} which uses Markov Chain Monte Carlo methods to fit ordinal probit regression models.

An overview on ordinal regression models in other statistical software packages like \proglang{Stata} \citep{Stata}, \proglang{SAS} \citep{SAS} or \proglang{SPSS} \citep{SPSS} is provided by \cite{Liu09}. These software packages include the \proglang{Stata} procedure \code{OLOGIT}, the \proglang{SAS} procedure \code{PROC LOGISTIC} and the \proglang{SPSS} procedure \code{PLUM} which perform ordinal logistic regression models. The software procedure \code{PLUM} additionally includes other link functions like probit complementary log-log, cauchit abd negative log-log. Ordinal models for multinomial data are available in the \proglang{SAS} package \code{PROC GENMOD}, while another implementation of ordinal logistic regression is available in \proglang{JMP} \citep{JMP}. In \proglang{Python} \citep{Python}, package \pkg{mord} \citep{mordpkg} implements ordinal regression methods.

While there are sufficient software tools in \proglang{R} which deal
with the univariate case, the ready-to-use packages for dealing with
the multivariate case fall behind, mainly due to computational
problems or lack of flexibility in the model specification. However,
there are some \proglang{R} packages which support correlated ordinal
data. One-dimensional normally distributed random effects in ordinal
regression can be handled by the \code{clmm()} function of package
\pkg{ordinal} \citep{ordinalR}; multiple possibly correlated random
effects are implemented in package \pkg{mixor} \citep{mixor}. Note
that this package uses multidimensional quadrature methods and
estimation becomes infeasible for increasing dimension of the random
effects. Bayesian multilevel models for ordinal data are implemented
in the package \pkg{brms} \citep{brmspkg}. Multivariate ordered
probit models, where the subject errors are assumed to follow a
multivariate normal distribution with a general correlation matrix,
can be estimated with package \pkg{PLordprob} \citep{PLordprob}, which
uses maximum composite likelihood methods estimation. This package
works well for standard applications but lacks flexibility. For
example, the number of levels of the ordinal responses needs to be
equal across all dimensions, threshold and regression coefficients are
the same for all multiple measurements and it does not account for
missing observations in the outcome variable.  None of these packages
support at the time of writing covariate dependent error structures. A
package which allows for different error structures in non-linear
mixed effects models is package \pkg{nlme} \citep{nlme}, even
though models dealing with ordinal data are not supported.

% Motivation
% ----------
The motivation of this package lies in a credit risk application,
where multiple credit ratings are assigned by various credit rating
agencies (CRAs) to firms over several years. CRAs have an important
role in financial markets, as they deliver subjective assessments or
opinions of an entity's creditworthiness, which are then used by the other players on the
market, such as investors and regulators, in their decision making
process.  Entities are assigned to rating classes by CRAs on an
ordinal scale by using both quantitative and qualitative criteria.
Ordinal credit ratings can be seen as a coarser version of an
underlying continuous latent process, which measures
creditworthiness. In the literature, this latent variable motivation
has been used in various credit rating models \citep[e.g.,][]{Blume98,
  afonso2009, Alp13, reusens2017sovereign}.

This setting is an example of an application where correlated ordinal
data arises naturally. On the one hand, multiple ratings for one firm
at the same point in time can be assumed to be correlated and on the
other hand, given the longitudinal dimension of the data, for each
rater, there is serial dependence in the ratings assigned over several
periods.  Moreover, aside from the need of a model class that can
handle correlated ordinal data, additional flexibility is desired due
to the following characteristics of the problem at hand:
% Firstly
Firstly, there is heterogeneity in the rating methodology. Raters use
different labeling as well as a different number of rating classes.
% Secondly
Secondly, the credit risk measure employed in assessing
creditworthiness can differ among raters (e.g., probability of default
versus recovery in case of default), which leads to heterogeneity in
the covariates, as raters might use different variables in their
rating process and assign different importance to the variables
employed.
% Thirdly
Thirdly, the data has missing values and is unbalanced, as firms can
leave the data set before the end of the observation period due to
various reasons such as default but also because of mergers and
acquisitions, privatizations, etc., or ratings can be withdrawn.
Moreover, there are missings in the multiple ratings, as not all firms
are rated by all CRAs at each time point.

% Contribution
% ------------
The paper and package \pkg{mvord} \citep{mvordpkg} for
\proglang{R} aim at providing a flexible framework for analyzing
correlated ordinal data by means of the class of multivariate ordinal
regression models. We offer the following features which (to the best
of our knowledge) enhance the currently available software for
multivariate ordinal regression models in \proglang{R}:
\begin{enumerate}[i]
\item Different error structures are available (such as a general, an
  autoregressive order one, or an equicorrelation structure for the
  subject error);
\item we also account for heterogeneity in the error structure among
  the subjects by allowing the use of subject-specific covariates in
  the specification of the error structure;
\item we implement a multivariate logit link for the class of
  multivariate ordinal regression models;
\item we allow for outcome specific threshold parameters;
\item we allow for outcome specific regression parameters;
\item the user can impose further restrictions on the threshold and
  regression parameters in order to achieve a more parsimonious model
  (e.g., using one set of thresholds for all outcome dimensions);
\item we offer the possibility to switch between different
  parameterizations, which are needed in ordinal models to ensure
  identifiability.
\end{enumerate}
%% Structure of the paper
This paper is organized as follows: Section~\ref{sect:model} provides
an overview of the model class and the estimation procedure, including
model specification and identifiability issues.
Section~\ref{sect:implementation} presents the main functions of the
package.  A description of the illustrative data sets used for
exemplifying the functionalities of the package is given in
Section~\ref{sect:data}. A couple of worked examples are given in
Section~\ref{sect:examples}. Section~\ref{sect:conclusion} concludes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model class and estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model class and estimation}\label{sect:model}
Multivariate ordinal regression models are an appropriate modeling
choice when a vector of correlated ordinal response variables,
together with covariates, is observed for each unit or subject in the
sample.  The response vector can be composed of different variables,
i.e., multiple measurements on the same subject (e.g., different
credit ratings assigned to a firm by different CRAs,
different survey questions answered by an interviewee, etc.) or
repeated measurements on the same variable at different time points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Model specification}\label{subsect:model}
Let $Y_{ij}$ denote the ordinal observation and $\bm x_{ij}$ be a
$p$ dimensional vector of covariates for subject~$i$ and outcome~$j$,
where $i=1,\, \dots,\, n$ and $j\in J_i$, for $J_i$ a subset of all
available outcomes $J$ in the data set. Moreover, we denote by $q=|J|$
and $q_i = |J_i|$ the number of elements in the set $J$ and $J_i$,
respectively. Following the cumulative link modeling approach
\citep{Mccullagh80}, the ordinal response $Y_{ij}$ is assumed to be a
coarser version of a latent continuous variable
$\widetilde Y_{ij}$.  The observable categorical outcome $Y_{ij}$ and
the unobservable latent variable $\widetilde Y_{ij}$ are connected by:
\begin{align*}
	Y_{ij} = r_{ij} \quad \Leftrightarrow \quad \theta_{j,r_{ij}-1} <
        \widetilde{Y}_{ij} \leq \theta_{j, r_{ij}}, \qquad r_{ij} \in
        \{1,\, \dots,\, K_j\},
\end{align*}
where $r_{ij}$ is a category out of $K_j$ ordered categories and $\bm
\theta_{j}$ is a vector of suitable threshold parameters for
outcome~$j$ with the following restriction: $ - \infty < \theta_{j,1}<
\dots < \theta_{j,K_j-1} < \infty$. Note that in this setting binary
observations can be treated as ordinal observations with two categories
($K_j = 2$).

The following linear model is assumed for the relationship between the
latent variable $\widetilde{Y}_{ij}$ and the vector of covariates
$\bm{x}_{ij}$:
\begin{align}\label{eqn:model}
   	\widetilde{Y}_{ij} &= \beta_{j0} + \bm{x}_{ij}^\top
        \bm{\beta}_j + \epsilon_{ij},
\end{align}
where $\beta_{j0}$ is an intercept term, $\bm{\beta}_j = (\beta_{j1},\,
\dots,\, \beta_{jp})^\top$ is a vector of regression coefficients, both
corresponding to outcome~$j$, and $\epsilon_{ij}$ is a mean zero error
term. We further assume the $n$~subjects to be independent and that the error
terms are uncorrelated with the covariates. Note that the number of ordered categories $K_j$ as well as the
threshold parameters $\bm\theta_j$ and the regression coefficients
$\bm\beta_j$ are allowed to vary across outcome dimensions $j \in J$
to account for possible heterogeneity across the response variables.

% Partial Prop. Odds
\paragraph{Category-specific regression coefficients}
By employing one set of regression coefficients $\bm\beta_j$
for all categories of the $j$-th outcome it is implied that the relationship
between the covariates and the responses does not depend on the category.
This assumption is called parallel regression or proportional odds
assumption \citep{Mccullagh80} and can be relaxed for one or more covariates by
allowing the corresponding regression coefficients to be category specific \citep[see e.g.,][]{Peterson90}.

%% Links
\paragraph{Link functions}
The dependence among the different responses is accounted for by
assuming that, for each subject~$i$, the vector of error terms
$\bm{\epsilon}_i = [\epsilon_{ij}]_{j\in J_i}$ follows a suitable
multivariate distribution. We consider two multivariate distributions
which correspond to the multivariate probit and logit link
functions. For the \emph{multivariate probit link} we assume that the
errors follow a multivariate normal distribution: $\bm{\epsilon}_i
\sim \mathcal{N}_{q_i}(\bm 0,\bm\Sigma_i)$. A \emph{multivariate
  logit link} is constructed by employing a multivariate logistic
distribution family with univariate logistic margins and a $t$ copula
with certain degrees of freedom.  For a vector $\bm z=(z_1,\,\dots,\,
z_q)^\top$, the multivariate logistic distribution function with $\nu$
degrees of freedom, mean $\bm \mu$ and covariance matrix $\bm \Sigma$
is defined as:
\begin{align}\label{eqn:logistic}
  F_{\nu, \bm\mu,\bm\Sigma}(\bm z) = t_{\nu, \bm
    R}(\{g_\nu((z_1-\mu_1)/\sigma_1),\, \dots,\,
  g_\nu((z_q-\mu_q)/\sigma_q)\}^\top),
\end{align}
where $t_{\nu,\bm R}$ is the $q$~dimensional multivariate
$t$~distribution with $\nu$ degrees of freedom and correlation matrix
$\bm R$ corresponding to $\bm \Sigma$, $g_\nu(x) =
t^{-1}_{\nu}(\exp(x)/(\exp(x) + 1))$, $t^{-1}_{\nu}$ is the quantile
function of the univariate $t$ distribution with $\nu$ degrees of
freedom and $\sigma^2_1,\,\dots,\,\sigma^2_q$ are the diagonal
elements of $\bm\Sigma$.

This $t$~copula based multivariate logistic family was proposed by
\cite{o2004bayesian} and later also employed by \cite{BIOM:BIOM12414}
and \cite{epubwu5389}, who approximate the distribution by a
multivariate $t$ distribution with the degrees of freedom chosen
appropriately. The employed distribution family differs from the
conventional multivariate logistic distributions of
\cite{gumbel1961bivariate} or \cite{malik1973multivariate} in that it
offers a more flexible dependence structure through the correlation
matrix of the $t$~copula, while still keeping the log odds
interpretation of the regression coefficients through the univariate
logistic margins.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IDENTIFIABILTY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Identifiability issues}\label{subsect:ident}
As the absolute scale and the absolute location are not identifiable
in ordinal models, further restrictions on the parameter set need to be
imposed.  Assuming $\bm\Sigma_i$ to be a covariance matrix with diagonal
elements $[\sigma^2_{ij}]_{j \in J_i}$, only the quantities
$\bm\beta_j/\sigma_{ij}$ and $(\theta_{j,r_{ij}}
-\beta_{j0})/\sigma_{ij}$ are identifiable in the model in
Equation~\ref{eqn:model}. Hence, in order to obtain an identifiable model the parameter set is typically
constrained in one of the following ways:
\begin{itemize}
  \item Fixing the intercept $\beta_{j0}$ (e.g., to zero), using
    flexible thresholds $\bm\theta_{j}$ and fixing $\sigma_{ij}$ (e.g.,
    to unity)  $\forall j \in J_i$, $\forall i \in \{1,\, \dots,\, n\}$;
  \item Leaving the intercept $\beta_{j0}$ unrestricted, fixing one
    threshold parameter (e.g., $\theta_{j,1}=0$) and fixing
    $\sigma_{ij}$ (e.g., to unity) $\forall j \in J_i$, $\forall i \in \{1,\, \dots,\, n\}$;
   \item Fixing the intercept $\beta_{j0}$ (e.g., to zero), fixing one
     threshold parameter (e.g., $\theta_{j,1}=0$) and leaving
     $\sigma_{ij}$ unrestricted $\forall j \in J_i$, $\forall i \in \{1,\, \dots,\, n\}$;
   \item Leaving the intercept $\beta_{j0}$ unrestricted, fixing two
     threshold parameters (e.g., $\theta_{j,1}=0$ and
     $\theta_{j,2}=1$) and leaving $\sigma_{ij}$ unrestricted $\forall
     j \in J_i$, $\forall i \in \{1,\, \dots,\, n\}$.\footnote{Note that this
     parameterization cannot be applied to the binary case.}
\end{itemize}
Note that the first two options are the most commonly used in the
literature. All of these alternative model parameterizations are
supported by the \pkg{mvord} package, allowing the user to choose the
most convenient one for each specific
application. Table~\ref{tab:Identifiability} in Section~\ref{sec:constraints.thresholds} gives an overview on
the identifiable parameterizations implemented in the package.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Error structures}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Basic model}
 The \emph{basic} multivariate ordinal regression model assumes that
 the correlation (and possibly variance, depending on the
 parameterization) parameters in the distribution function of the
 $\bm\epsilon_i$'s are constant for all subjects~$i$.%=1,\dots,n$.
 %% Correlation
 \paragraph{Correlation}
 The dependence between the multiple measurements or outcomes can be
 captured by different correlation structures. Among them, we
 concentrate on the following three:
 \begin{itemize}
   \item the \emph{general} correlation structure assumes different correlation
     parameters between pairs of outcomes $corr(\epsilon_{ik},
     \epsilon_{il}) = \rho_{kl}$;
   \item the \emph{equicorrelation} structure $corr(\epsilon_{ik},
     \epsilon_{il}) = \rho$ implies that the correlation between
     all pairs of outcomes is constant;
   \item when faced with longitudinal data, especially when moderate
     to long subject-specific time series are available, an \emph{$AR(1)$
     autoregressive} correlation model of order one can be employed. For given equally spaced time points $t_1$, $t_2$, $\dots$, $t_q$
     this $AR(1)$ error structure implies an exponential decay in the correlation with the
     lag. If $t_{k}$ and $t_{l}$ are the time points when $Y_{ik}$
     and $Y_{il}$ are observed, then $corr(\epsilon_{ik},
     \epsilon_{il}) = \rho^{|t_{k}-t_{l}|}$.
 \end{itemize}
 %% Variance
 \paragraph{Variance}
  If a parameterization which supports the estimation of the variance
  of the latent processes is used (see Section~\ref{subsect:ident}),
  it is assumed that $\VAR(\epsilon_{ij}) = \sigma^2_{j}$.

 \subsubsection{Extending the basic model}
 In some applications, the constant correlation (and variance)
 structure across subjects may be too restrictive. We hence extend the
 basic model by allowing the use of covariates in the
 correlation (and variance) specifications.

\paragraph{Correlation}
  For each subject~$i$ and each pair $(k,l)$ from the set $J_i$, the
  correlation parameter $\rho_{ikl}$ is assumed to depend on a vector
  $\bm{s}_{i}$ of $m$ subject-specific covariates. The hyperbolic
  tangent transformation allows us to reparameterize the linear term
  $\alpha_{0kl} + \bm{s}_{i}^\top \bm \alpha_{kl}$ in terms of a
  correlation parameter:
  \begin{align*}
  \frac{1}{2} \log \left( \frac{1+ \rho_{ikl}}{1-\rho_{ikl}}\right) =
  \alpha_{0kl} + \bm{s}_{i}^\top \bm\alpha_{kl}, \qquad \rho_{ikl} =
  \frac{e^{2 (\alpha_{0kl} + \bm{s}_{i}^\top \bm\alpha_{kl})} -
    1}{e^{2 (\alpha_{0kl} + \bm{s}_{i}^\top \bm\alpha_{kl})} + 1}.
  \end{align*}
 If $\bm{\alpha}_{kl}=0$ for all $k,l\in J_i$, this model would
 correspond to the general correlation structure in the basic
 model. Moreover, if $\alpha_{0kl}=0$ and $\bm{\alpha}_{kl}=0$ for all $k,l\in
 J_i$, the correlation matrix is the identity matrix and the responses
 are uncorrelated.

For the more parsimonious error structures of equicorrelation and
$AR(1)$, in the extended model the correlation parameters are modeled
as:
\begin{align*}
  \frac{1}{2} \log \left( \frac{1+ \rho_{i}}{1-\rho_{i}}\right) =
  \alpha_{0} + \bm{s}_{i}^\top \bm{\alpha}, \qquad \rho_{i} =
  \frac{e^{2 (\alpha_{0} + \bm{s}_{i}^\top \bm{\alpha})} - 1}{e^{2
      (\alpha_{0} + \bm{s}_{i}^\top \bm{\alpha})} + 1}.
\end{align*}

\paragraph{Variance}
Similarly, one could model the heterogeneity among the subjects through the variance parameters $\VAR(\epsilon_{ij}) = \sigma^2_{ij}$
by employing the following linear model on the log-variance:
\begin{align*}
  \text{log}(\sigma^2_{ij}) = \gamma_{0j} + \bm{s}_{i}^\top \bm\gamma_{j}.
\end{align*}

The correlation (or covariance) matrix $\bm\Sigma_i$ must be
positive-semi-definite. This can be ensured by the use of special
algorithms such as the one proposed by \cite{higham1988computing}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Composite Likelihood Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Composite likelihood estimation} \label{sect:estimation}
In order to estimate the model parameters we use a composite
likelihood approach, where the full likelihood is approximated by a
pseudo-likelihood which is constructed from lower dimensional
marginal distributions, more specifically by ``aggregating'' the
likelihoods corresponding to pairs of observations \citep{varin_overview}.

For a given parameter vector $\bm\delta$, which contains the threshold
parameters, the regression coefficients and the parameters of the
error structure, the likelihood is given by:
\begin{align*}
\mathscr{L} (\bm\delta) &= \prod_{i=1}^n \Prob\bigg(\bigcap_{j\in
  J_{i}}\{Y_{ij} = r_{ij}\}\bigg)^{w_{i}} =\prod_{i=1}^n \bigg(
\int_{D_{i}} f_{q_{i}}(\bm{\widetilde{Y}}_{i}; \bm
\delta)d^{q_{i}}\widetilde{\bm{Y}}_{i}\bigg)^{w_{i}},
 \end{align*}
where $D_{i} =
\prod_{j\in J_{i}} (\theta_{j,r_{ij}-1}, \theta_{j,r_{ij}}]$ is a
  Cartesian product, $w_{i}$ are subject-specific non-negative
  weights, which are set to one in the default case, and $f_{q_{i}}$
  is the $q_{i}$-dimensional density of the error terms
  $\bm\epsilon_i$. We approximate this full likelihood by a pairwise likelihood which is
constructed from bivariate marginal distributions. If the number of
observed outcomes for subject~$i$ is less than two ($q_i<2$), the
univariate marginal distribution enters the likelihood. The pairwise
log-likelihood function is obtained by:
\begin{align}\label{eqn:logpl}
  p\ell(\bm\delta)= \sum_{i=1}^n w_i \biggl[& \mathbbm{1}_{\{q_i \geq
      2\}}\sum_{\substack{k<l\\ k,l\in J_i}}\log\left(\Prob(Y_{ik} = r_{ik},
    Y_{il} = r_{il})\right)+\nonumber\\ & \mathbbm{1}_{\{q_i = 1\}}
    \mathbbm{1}_{\{k \in J_i\}}\biggl.\log\left(\Prob(Y_{ik} = r_{ik})
    \right)\biggr].
\end{align}
Denoting by $U_{ij} = ( \theta_{j,r_{ij}} - \beta_{j0} -\bm
x_{ij}^\top\bm{\beta}_j)/\sigma_{ij}$ the upper and by $L_{ij} =
(\theta_{j,r_{ij}-1} - \beta_{j0} -
\bm{x}_{ij}^\top\bm{\beta}_j)/\sigma_{ij}$ the lower integration
bounds and by $f_1$ and $f_2$ the uni- and bivariate density functions
corresponding to the error distribution, the uni- and bivariate
probabilities are given by:
\begin{align*}
\Prob(Y_{ik} &= r_{ik}, Y_{il} = r_{il})
=\displaystyle\int_{L_{ik}}^{U_{ik}}
\displaystyle\int_{L_{il}}^{U_{il}} f_{2}(v_{ik},v_{il}; \rho_{ikl})
dv_{ik}dv_{il},\\ \Prob(Y_{ik} &= r_{ik}) =\displaystyle
\int_{L_{ik}}^{U_{ik}} f_{1}(v_{ik}) dv_{ik}.
\end{align*}

The maximum pairwise likelihood estimates
$\hat{\bm\delta}_{\text{PL}}$ are obtained by direct maximization of
the composite likelihood given in Equation~\ref{eqn:logpl}. The
threshold and error structure parameters to be estimated are
reparameterized such that unconstrained optimization can be
performed. Firstly, we reparameterize the threshold parameters in
order to achieve monotonicity. Secondly, for all unrestricted
correlation (and covariance) matrices we use the spherical
parameterization of \cite{Pinheiro96}. This parameterization has the
advantage that it can be easily applied to correlation
matrices. Thirdly, if we assume to have equicorrelated or $AR(1)$
errors, we use the hyperbolic tangent transformation.

Computation of the standard errors is needed in order to quantify the
uncertainty of the maximum pairwise likelihood estimates. Under
certain regularity conditions, the maximum pairwise likelihood
estimates are consistent as the number of responses is fixed and
$n\rightarrow \infty$. In addition, the maximum pairwise likelihood
estimator is asymptotically normal with asymptotic mean $\bm\delta$
and a covariance matrix which equals the inverse of the Godambe
information matrix:
\begin{align*}
  G(\bm\delta)^{-1} = H^{-1}(\bm\delta)V(\bm\delta) H^{-1}(\bm\delta),
\end{align*}
where $H(\bm\delta)$ the Hessian (sensitivity matrix) and $V(\bm\delta)$ the
variability matrix. The variability matrix $V(\bm\delta)$ and the
Hessian $H(\bm\delta)$ can be estimated as follows:
\begin{align*}
\hat V(\bm\delta)= \frac{1}{n} \sum_{i=1}^n \left(\frac{\partial
  p\ell_i(\hat{\bm\delta}_{\text{PL}})}{\partial \bm\delta}\right)
\left(\frac{\partial p\ell_i(\hat{\bm\delta}_{\text{PL}})}{\partial
  \bm\delta}\right)^{\top},
\end{align*}
\begin{align*}
\hat H(\bm\delta) &= -\frac{1}{n} \sum_{i=1}^n \frac{\partial^2
  p\ell_i(\hat{\bm\delta}_{\text{PL}})}{\partial
  \bm\delta\partial\bm\delta^\top}\\ &=\frac{1}{n} \sum_{i=1}^n
\sum_{\substack{k<l\\ k,l\in J_i}} \left(\frac{\partial
  p\ell_{ikl}(\hat{\bm\delta}_{\text{PL}})}{\partial \bm\delta}\right)
\left(\frac{\partial
  p\ell_{ikl}(\hat{\bm\delta}_{\text{PL}})}{\partial
  \bm\delta}\right)^\top,
\end{align*}
where $p\ell_i(\bm\delta)$ is the component of the pairwise
log-likelihood corresponding to subject~$i$ and $p\ell_{ikl}(\bm\delta)$ corresponds to subject~$i$ and pair $(k,l)$.

In order to compare different models, the composite likelihood
information criterion by \cite{varin2005note} can be used:
$\text{CLIC}(\bm\delta) = -2\ p\ell(\hat{\bm\delta}_{\text{PL}}) +
k\ \mathrm{tr}(\hat V(\bm\delta)\hat H(\bm\delta)^{-1})$ (where $k=2$
corresponds to CLAIC and $k=\text{log}(n)$ corresponds to CLBIC). A
comprehensive overview and further details on the properties of the
maximum composite likelihood estimates is provided in \cite{Varin08}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interpretation of the coefficients}
Unlike in linear regression models, the interpretation of the
regression coefficients and of the threshold parameters in ordinal
models is not straightforward. Estimated thresholds and coefficients
represent only signal to noise ratios and cannot be interpreted
directly (see Section~\ref{subsect:ident}). For one particular
outcome~$j$, the coefficients can be interpreted in the same way as in
univariate cumulative link models.  Let us assume without loss of
generality that a higher latent score leads to better ratings on the
ordinal scale. This implies that the first category is the worst and
category $K_j$ is the best category.  The marginal cumulative
probabilities implied by the model in Equation~\ref{eqn:model} are
then given by the following relationship:
\begin{align*}
\Prob(Y_{ij} \leq r_{ij}|\bm x_{ij}) = \Prob(\bm{x}_{ij}^\top \bm{\beta}_j +
\epsilon_{ij} \leq \theta_{j,r_{ij}}) = \Prob(\epsilon_{ij} \leq \theta_{j,r_{ij}} -
\bm{x}_{ij}^\top \bm{\beta}_j)  = F_1(\theta_{j,r_{ij}} - \bm{x}_{ij}^\top
\bm{\beta}_j),
\end{align*}
where $F_1$ is the marginal (univariate) distribution function of the
errors for the $j$-th outcome.

%marginal/partial effects
One natural way to interpret ordinal regression models is to analyze
partial effects, where one is interested in how a marginal change in one
variable $x_{ijv}$ changes the outcome distribution. The partial
probability effects in the cumulative model are given by:
\begin{align*}
  \delta_{r,v}^j(\bm x_{ij}) = \frac{\partial \Prob(Y_{ij} = r_{ij}|
    \bm x_{ij})}{\partial x_{ijv}} = -\left(f_1(\theta_{j,r_{ij}} -
  \bm{x}_{ij}^\top \bm{\beta}_j) - f_1(\theta_{j,r_{ij}-1} -
  \bm{x}_{ij}^\top \bm{\beta}_j)\right)\beta_{jv},
\end{align*}
where $f_1$ is the density corresponding to $F_1$, $x_{ijv}$ is the
$v$-th element in $\bm x_{ij}$ and $\beta_{jv}$ is the $v$-th element
in $\bm \beta_{j}$. In case of discrete variables it is more
appropriate to consider the changes in probability before and after
the change in the variable instead of the partial effects using:
\begin{align*}
\Delta \Prob(Y_{ij} = r_{ij}|\bm x_{ij}, \tilde{\bm x}_{ij}) =
\Prob(Y_{ij} = r_{ij}|\tilde{\bm x}_{ij}) - \Prob(Y_{ij} = r_{ij} |
\bm x_{ij}),
\end{align*}
where all elements of $\tilde{\bm x}_{ij}$ are equal to $\bm x_{ij}$
except for the $v$-th element, which is equal to $\tilde{x}_{ijv} =
x_{ijv} + \Delta x_{ijv}$ for the discrete change $\Delta x_{ijv}$ in
the variable $x_v$. We refer to \cite{Hensher10} and \cite{Boes06} for
further discussion of the interpretation of partial effects in ordered
response models.

In the presence of the probit link function, we have the following
relationship between the cumulative probabilities and the latent
process:
\begin{align*}
\Phi^{-1}\left(\Prob(Y_{ij} \leq r_{ij}|\bm{x}_{ij})\right) = \theta_{j,r_{ij}} -
\bm{x}_{ij}^\top \bm{\beta}_j.
\end{align*}

An increase of one unit in a variable $x_{jv}$ (given that all other
variables are held constant) changes the probit of the probability
that category $r$ or lower is observed by the value of the coefficient
$\beta_{jv}$ of this variable. In other words $\Prob(Y_{ij} \leq
r_{ij}|\bm{x}_{ij})$, the probability that category $r_{ij}$ or lower
is observed, changes by the increase/decrease in the distribution
function.  Moreover, predicted probabilities for all ordered response
categories can be calculated and compared for given sets of
explanatory variables.

In the presence of the logit link function, the regression
coefficients of the underlying latent process are scaled in terms of
marginal log odds \citep{Mccullagh80}:
\begin{align*}
  \log\left(\frac{\Prob(Y_{ij} \leq r_{ij}|\bm{x}_{ij})}{\Prob(Y_{ij}
    > r_{ij}|\bm{x}_{ij})} \right) = \theta_{j,r_{ij}} -
  \bm{x}_{ij}^\top \bm{\beta}_j.
\end{align*}
For a one unit increase in one variable $x_{jv}$ holding all the
others constant, we expect a change of size of the coefficient
$\beta_{jv}$ of this variable in the expected value on the log odds
scale.  Due to the fact that the marginal effects of the odds ratios
do not depend on the category, one often exponentiates the
coefficients in order to obtain the following convenient
interpretation in terms of odds ratios:
\begin{align*}
\frac{\Prob(Y_{ij} \leq r_{ij}|\bm{x}_{ij})/\Prob(Y_{ij} >
r_{ij}|\bm{x}_{ij})}{\Prob(Y_{ij} \leq r_{ij}|\tilde{\bm{x}}_{ij})/\Prob(Y_{ij}
> r_{ij}|\tilde{\bm{x}}_{ij})} = \exp((\tilde{\bm{x}}_{ij} -
\bm{x}_{ij})^\top \bm{\beta}_j).
\end{align*}
This means for a one
unit increase in $x_{jv}$ holding all the other variables constant, changes the
odds ratio by $e^{\beta_{jv}}$. In other words, the odds after a one unit
change in $x_{jv}$ are the odds before the change multiplied by $e^{-\beta_{jv}}$
\begin{align*}
  \frac{\Prob(Y_{ij} \leq r_{ij}|\bm{x}_{ij})}{\Prob(Y_{ij} > r_{ij}|\bm{x}_{ij})}
\exp(-\bm{\beta}_j) = \frac{\Prob(Y_{ij} \leq
r_{ij}|\tilde{\bm{x}}_{ij})}{\Prob(Y_{ij} > r_{ij}|\tilde{\bm{x}}_{ij})}.
\end{align*}

If the regression coefficients vary across the multiple responses,
they cannot be compared directly due to the fact that the measurement
units of the underlying latent processes differ. Nevertheless, one
possibility to compare coefficients is through concept of
\textit{importance}. \cite{reusens2017sovereign} extend an approach
for comparing coefficients of probit and logit models by
\cite{Hoetker07} in order to compare the coefficients across repeated
measurements. They analyze the \textit{importance ratio}
\begin{align*}
  R_{jv} = \frac{\beta_{jv}}{\beta_{j,base}},
\end{align*}
where $\beta_{j,base}$ is the coefficient of a \textit{base} variable and $v$ is one of the remaining $p-1$ variables. This ratio can be interpreted as follows: A one unit increase in the variable $v$ has in expectation the same effect in the $base$ variable multiplied by the ratio $R_{jv}$.
Another interpretation is the so called \textit{compensation variation}: %\citep{Boes06, Train98, Train03}
The ratio is the required increase in the $base$ variable that is necessary to compensate a one unit decrease in the variable $v$ in a way that the score of the outcome remains the same.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IMPLEMENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}\label{sect:implementation}
Multivariate ordinal regression models in the \proglang{R} package
\pkg{mvord} are fitted using the function \code{mvord()}:
\begin{verbatim}
mvord(formula, data, error.structure = cor_general(~ 1), link = mvprobit(),
  index = NULL, response.names = NULL, response.levels = NULL,
  threshold.constraints = NULL, threshold.values = NULL,
  coef.constraints = NULL, coef.values = NULL, weights = NULL, se = TRUE,
  offset = NULL, scale = FALSE, start.values = NULL, solver = "newuoa",
  PL.lag = NULL, control = list(maxit=200000, trace = 1, kkt = FALSE))
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Data structure}
We use a long format for the input of \code{data}, where each row
contains a subject index~$i$ (\code{subject_index}), a multiple measurement
index~$j$ (\code{multiple_measurement_index}), an ordinal response (\code{Y}) and
all the covariates (\code{X1} to \code{Xp})\footnote{If the covariates have different scales, numerical
instabilities may occur in the estimation procedure. In such cases we suggest to e.g., standardize the covariates or transform them to have a similar scale.}. This long format data
structure is internally transformed to an $n\times q$ matrix of responses $\bm Y_1$ to $\bm Y_q$ (which
contains \code{NA} in the case of missing entries) and a list of covariate
matrices $\bm X_j$ for all $j \in J$. In order to construct these objects,
the subject index~$i$ and the multiple measurement index~$j$ should be
specified. This can be performed by an optional argument \code{index},
a character vector of length two, specifying the column names of the
subject index and the multiple measurement index in \code{data}.
\begin{verbatim}
index = c("subject_index", "multiple_measurement_index")
\end{verbatim}
The default value of \code{index} is \code{NULL} assuming that the
first column of \code{data} contains the subject index~$i$ and the second
column the multiple measurement index~$j$.

If specific constraints are imposed on the threshold parameters and/or
on the regression coefficients, it is important to know which level of
the multiple measurement index~$j$ corresponds to the first dimension,
second dimension and so on. Hence, a well defined index $j \in J$ for
the multiple measurements is needed. Therefore, a vector
\code{response.names} is used to define the index number of the
multiple measurement. If we assume to have $q=4$ outcomes \code{Y1} to \code{Y4}, the following \code{response.names} are set:
\begin{verbatim}
response.names = c("Y1", "Y2", "Y3", "Y4")
\end{verbatim}
The default value of \code{response.names} is \code{NULL} giving the natural ordering of the levels of the factor variable for all the multiple measurements. The ordering of \code{response.names} always specifies the index of the multiple measurement unit $j \in J$. This ordering is essential when putting constraints on the parameters and when setting \code{response.levels}.
% #
% # response.levels <-  list(Y1 = rev(LETTERS[1:6]),
% #                          Y2 = rev(LETTERS[1:6]),
% #                          Y3 = rev(LETTERS[7:13]),
% #                          Y4 = c("O", "N"))
% # response.levels
% # @
\begin{verbatim}
response.levels = list(Y1 = c("F", "E", "D", "C", "B", "A"),
                       Y2 = c("F", "E", "D", "C", "B", "A"),
                       Y3 = c("M", "L", "K", "J", "I", "H", "G"),
                       Y4 = c("O", "N"))
\end{verbatim}
If the categories differ across multiple measurements (either the number of categories and/or the category labels) one needs to specify the \code{response.levels} explicitly. This is performed by a list of length $q$, where each element contains the names of the levels of the ordered categories in ascending (or if desired descending) order.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\subsection{Formula}
The ordinal responses \code{Y} and the $p$~covariates \code{X1}, \code{X2}, $\dots$, \code{Xp}    are passed by a
\code{formula} object. Intercepts can be included or excluded in the
model depending on the model parameterization chosen in order to
ensure identifiability:
\paragraph{Model without intercept}
If the intercept should be removed, the \code{formula} has
the following form:
\begin{verbatim}
formula = Y ~ 0 + X1 + ... + Xp
\end{verbatim}
\paragraph{Model with intercept}
If one wants to include an intercept in the model, there are two
equivalent possibilities to set the model \code{formula}. Either the
intercept is included explicitly by:
\begin{verbatim}
formula = Y ~ 1 + X1 + ... + Xp
\end{verbatim}
or by
\begin{verbatim}
formula = Y ~ X1 + ... + Xp
\end{verbatim}
\paragraph{Note on intercept in formula}
We differ in our approach of specifying the model \code{formula}
from the formula objects in e.g., \code{MASS::polr()} or \code{ordinal::clm()}, in that we allow the user to specify models without intercept. This option is not supported in the \pkg{MASS} and \pkg{ordinal} packages, where an intercept is always specified in \code{formula} as the threshold parameters are treated as intercepts.
We choose to allow for this option, in order to have a correspondence
to the identifiability constraints presented in Section~\ref{subsect:ident}.

Even so, the user should be aware that
the threshold parameters are basically category and outcome specific intercepts.
This implies that, even if the intercept is explicitly removed from the model
through the \code{formula} object and hence set to zero, the rest of the
covariates should be specified in such a way that multicolinearity does not arise (as if an intercept were specified). This is of primary importance when including categorical covariates, where one category will be taken as baseline by default.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Link functions}
The multivariate link functions are specified as objects of class
\code{`mvlink'}.  An object of class \code{`mvlink'} is a list with
elements specifying the distribution function of the errors, functions
for computing the corresponding univariate and bivariate
probabilities, as well as additional arguments specific to each
link. If gradient functions passed on, these will be used in the
computation of the standard errors.  This design has been inspired by
the design of the \code{`family'} class in package \pkg{stats} and it
facilitates the addition of new link functions to the package.

We offer two different multivariate link functions, the
multivariate probit link and a multivariate logit link.  For the
multivariate probit link a multivariate normal distribution for the
errors is applied. The normal bivariate probabilities which enter the
pairwise log-likelihood are computed with package \pkg{pbivnorm}
\citep{pbivnorm}. The multivariate probit link can be specified by:
\begin{verbatim}
link = mvprobit()
\end{verbatim}
For the multivariate logit link a $t$~copula based multivariate
distribution with logistic margins is used (as explained in
Section~\ref{subsect:model}) and can be specified by:
\begin{verbatim}
link = mvlogit(df = 8L)
\end{verbatim}
The \code{mvlogit()} function has an optional integer valued argument
\code{df} which specifies the degrees of freedom to be used for the
$t$~copula.  The default value of the degrees of freedom parameter is
8. When choosing $\nu \approx 8$, the multivariate logistic
distribution in Equation~\ref{eqn:logistic} is well approximated by a
multivariate $t$ distribution \citep{o2004bayesian}. This is also the
value chosen by \cite{BIOM:BIOM12414} in their
analysis. We restrict the degrees of freedom to be integer valued
because the most efficient routines for computing bivariate
$t$~probabilities do not support non-integer degrees of freedom. We
use the \proglang{Fortran} code from Alan Genz \citep{Genz09} to
compute the bivariate $t$~probabilities.  As the degrees of freedom
parameter is integer valued, we do not estimate it in the optimization
procedure. If the optimal degrees of freedom are of interest, we leave
the task of choosing an appropriate grid of values of \code{df} to the
user, who should then estimate a separate model for each value in the
grid. The best model can be chosen by CLAIC or CLBIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Error structures}
Different error structures are implemented in \pkg{mvord} and can be specified
through the argument \code{error.structure}. The error structure objects are
of class \code{`error_struct'}. This approach is slightly different than in package \pkg{nlme},
where the error structure is defined by two classes: \code{`varFunc'} for the variance function
and \code{`corStruct'} for the correlation structure.
We also define the following subclasses for the error structures: \code{`cor_general'}
(similar to  \pkg{nlme}'s \code{`corSymm'}),  \code{`cor_equi'} (similar to \code{`corCompSymm'}),
\code{`cor_ar1'}  (similar to \code{`corAR1'}) and
\code{`cov_general'} (similar to \code{`corSymm'} with variance function \code{`varIdent'}).

The different error structures are chosen in the following way:
\subsubsection{Basic model}
\paragraph{Correlation}
For the basic model specification the following correlation structures are implemented in \pkg{mvord}:
\begin{itemize}
  \item \code{cor_general()} -- A general error structure, where the
  correlation matrix of the error terms is unrestricted and constant
  across all subjects: $corr(\epsilon_{ik},
  \epsilon_{il})=\rho_{kl}$. This error structure is among the most
  common in the literature \citep[e.g.,][]{Scott02, Bhat10, Pagui15}.
\begin{verbatim}
error.structure = cor_general(formula = ~ 1)
\end{verbatim}
  \item \code{cor_equi()} -- An equicorrelation structure with
    $corr(\epsilon_{ik}, \epsilon_{il})=\rho$ is used.
\begin{verbatim}
error.structure = cor_equi(formula = ~ 1)
\end{verbatim}
\item \code{cor_ar1()} -- An autoregressive error structure of order one with
  $corr(\epsilon_{ik}, \epsilon_{il}) = \rho^{|t_{k}-t_{l}|}$ is used.
\begin{verbatim}
error.structure = cor_ar1(formula = ~ 1)
\end{verbatim}
\end{itemize}
\paragraph{Variance}
A model with variance parameters $\VAR(\epsilon_{ij})=\sigma^2_j$
corresponding to each outcome, when the identifiability requirements
are fulfilled, can be specified in the following way:
\begin{itemize}
  \item the estimation of $\sigma^2_j$ is only implemented in
    combination with the general correlation structure. Hence, the
    unrestricted covariance matrix of the error terms can be estimated
    by:
\begin{verbatim}
error.structure = cov_general(formula = ~ 1)
\end{verbatim}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Extending the basic model}
\paragraph{Correlation}
\begin{itemize}
  \item For the heterogeneous general correlation structure, the
    current implementation only allows the use of one \code{factor}
    variable \code{f} as covariate. As previously mentioned, this
    factor variable should be subject-specific and hence should not
    vary across the multiple responses.  This implies that a
    correlation matrix will be estimated for each factor level.
    %% only up to maximum 30 levels are allowed.
\begin{verbatim}
error.structure = cor_general(formula = ~ f)
\end{verbatim}
\item Estimating an equicorrelation structure depending on
  subject-specific covariates:
\begin{verbatim}
error.structure = cor_equi(formula = ~ S1 + ... + Sm)
\end{verbatim}
\item Estimating an $AR(1)$ correlation structure depending on
  subject-specific covariates:
\begin{verbatim}
error.structure =  cor_ar1(formula = ~ S1 + ... + Sm)
\end{verbatim}
\end{itemize}
\paragraph{Variance}
\begin{itemize}
\item As in the basic model, the estimation of the heterogeneous
  variance parameters can be performed for the
  general correlation structure. A subject-specific factor
  \code{f} can be used as covariate in the log variance equation. This
  can be specified by:
\begin{verbatim}
error.structure = cov_general(formula = ~ f)
\end{verbatim}
In addition to the correlation matrices, which are estimated for
each factor level of \code{f}, a vector of dimension $q$ of variance
parameters will be estimated for each factor level.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Constraints on thresholds} \label{sec:constraints.thresholds}
The package supports constraints on the threshold
parameters. Firstly, the user can specify whether the threshold
parameters should be equal across some or all response dimensions.
Secondly, the values of some of the threshold parameters can be fixed.
This feature is important for the users who wish to further restrict
the parameter space of the thresholds or who wish to specify values
for the threshold parameters other than the default values used in the
package.  Note that fixing some of the thresholds is needed for some
of the parameterizations presented in
Table~\ref{tab:Identifiability} in order to ensure identifiability
of the model.
%%%%%%%%%%%%%%%%%%
\subsubsection{Threshold constraints across responses}
Such constraints can be imposed by a vector of positive integers
\code{threshold.constraints}, where dimensions with equal threshold
parameters get the same integer. When restricting two outcome
dimensions to be equal, one has to be careful that the number of
categories in the two outcome dimensions must be the same. In an example with $q=4$ different outcomes, if one wishes to restrict the
threshold parameters of the first and second outcomes \code{Y1} and \code{Y2} to be equal ($\bm\theta_1 = \bm\theta_2$). These restrictions can be specified as:
% # <<
% # threshold.constraints <- c(1, 1, 2, 3)
% # names(threshold.constraints) <- paste0("Y", 1:4)
% # threshold.constraints
% # @
\begin{verbatim}
threshold.constraints = c(1, 1, 2, 3)
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Fixing threshold values}
Values for the threshold parameters can be specified by the argument
\code{threshold.values}. For this purpose the user can pass a
\code{list} with $q$ elements, where each element is a \code{vector}
of length $K_j - 1$ (where $K_j$ is number of ordered categories for
ordinal outcome~$j$). A numeric value in this vector fixes the
corresponding threshold parameter to the specified value while
\code{NA} leaves the parameter flexible and indicates it should be
estimated.

After specifying the error structure (through the
\code{error.structure} argument) and whether an intercept should be
estimated or not (in the \code{formula} argument), the user can choose
among five possible options for fixing the thresholds:
\begin{itemize}
  \item leaving all thresholds flexible;
  \item fixing the first threshold $\theta_{j,1}$ to a constant $a_j$ for all $j\in J$;
  \item fixing the first and second thresholds
    $\theta_{j,1}=a_j$, $\theta_{j,2}=b_j$ for all outcomes with $K_j > 2$;
  \item fixing the first and last thresholds
     $\theta_{j,1}=a_j$, $\theta_{j,K_j-1}=b_j$ for all outcomes with $K_j > 2$;
  \item an extra option is fixing all of the threshold parameters, for all $j\in J$.
\end{itemize}
Note that the option chosen needs to be consistent across the
different outcomes (e.g., it is not allowed to fix the first and the last
threshold for one outcome and the first and the second threshold for a
different). Table~\ref{tab:Identifiability} provides information
about the options available for each combination error structure and
intercept, as well as about the default values in case the user does
not specify any threshold values.

\begin{table}[h!]
\setlength{\tabcolsep}{2.5pt}
\begin{tabularx}{\textwidth}{|C|C|C|C|C|C|C|}%{p{1.4cm}|p{1.4cm}|p{1.3cm}|X|X|X|X|}%{|C|C|C|C|C|C|C|C|}
\hline
\multirow{4}{=}{\centering Error structure}&\multirow{4}{=}{\centering Intercept} & \multicolumn{5}{c|}{Threshold parameters}\\
\cline{3-7}
 & & all flexible & one fixed & two fixed & two fixed & all fixed \\
&   & & $\theta_{j,1}=a_j$ & $\theta_{j,1}=a_j$  &$\theta_{j,1}=a_j$  & \\
 &  &&& $\theta_{j,2}=b_j$ & $\theta_{j,K_j-1}=b_j$ & \\
\hline
\multirow{2}{*}{\code{cor}}  &no     &\textcolor{Green}{\checkmark} &\checkmark&\checkmark&\checkmark&\checkmark\\
 &yes    &   &\textcolor{Green}{\checkmark} & \checkmark &\checkmark&\checkmark\\
\hline
\multirow{2}{*}{\code{cov}}  & no &   & \textcolor{Green}{\checkmark } & \checkmark &  \checkmark &  \checkmark\\
 &yes&  &  & \textcolor{Green}{\checkmark }&\checkmark &\checkmark\\
\hline
\end{tabularx}
\caption{This table displays different model parameterizations in the presence of ordinal observations ($K_j > 2 \; \forall j \in J$). The
  row \code{cor} includes error structures \code{cor\_general},
  \code{cor\_equi} and \code{cor\_ar1}, while row \code{cov} includes the
  error structure \code{cov\_general}. The minimal restrictions (default) to
  ensure identifiability are given in green. The default threshold values (in
  case \code{threshold.values = NULL}) are always $a_j = 0$ and $b_j = 1$.}
  \label{tab:Identifiability}
\end{table}

In the presence of binary observations ($K_j = 2$), if a \code{cov_general} error structure is used,
the intercept has always to be fixed to some value due to identifiability constraints. In a correlation
structure setting no further restrictions are required.

For example, the following restrictions on the threshold parameters:
\begin{itemize}
\item $\theta_{11} = -4 \leq \theta_{12} \leq \theta_{13}\leq
  \theta_{14}\leq \theta_{15} \leq \theta_{16} $;
\item $\theta_{21} = -4 \leq \theta_{22} \leq \theta_{23}\leq
  \theta_{24}\leq \theta_{25} \leq \theta_{26} $;
\item $\theta_{31} = -5 \leq \theta_{32} \leq \theta_{33}\leq
  \theta_{34}\leq \theta_{35} \leq \theta_{36} \leq \theta_{37}$;
\item $\theta_{41} = 0$,
\end{itemize}
are implemented as:
% # <<
% # threshold.values <- list(c(-4, NA, NA, NA, NA, NA),
% #                          c(-4, NA, NA, NA, NA, NA),
% #                          c(-5, NA, NA, NA, NA, NA, NA),
% #                          c(0))
% # names(threshold.values) <- paste0("Y", 1:4)
% # threshold.values
% # @
\begin{verbatim}
threshold.values = list(c(-4, NA, NA, NA, NA, NA), c(-4, NA, NA, NA, NA, NA),
  c(-5, NA, NA, NA, NA, NA, NA), 0)
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Constraints on coefficients}
The package supports
constraints on the regression coefficients. Firstly, the
user can specify whether the regression coefficients should be equal
across some or all response dimensions. Secondly, the values of some
of the regression coefficients can be fixed.

As there is no unanimous way to specify such constraints, we offer
two options. The first option is similar to the specification of constraints on the thresholds.
The constraints can be specified in this case as a vector or matrix of integers, where coefficients getting same integer value are set equal.
Values of the regression coefficients can be fixed through a matrix.
Alternatively, constraints on the regression coefficients can be specified
by using the design employed by the \pkg{VGAM} package.
The constraints in this setting are set through a named list,
where each element of the list contains a matrix of full-column rank.
If the values of some regression coefficients should be fixed, offsets can be used.
This design has the advantage that it supports
constraints on outcome-specific as well as category-specific
regression coefficients. While the first option has the advantage of requiring a more concise input, it does not support category-specific coefficients.
The second option offers a more flexible design in this respect.

%
\subsubsection{Coefficient constraints across responses}
Such constraints can be specified by the argument
\code{coef.constraints}, which can be either a vector or a matrix of
integer values. If vector constraints of the type $\bm \beta_{k} = \bm \beta_{l}$, are desired, which
should hold for all $p$~regression coefficients corresponding to
outcome $k$ and $l$, the easiest way to specify this is by means of a
vector of integers of dimension~$q$, where outcomes with equal vectors
of regression coefficients get the same integer. For example, for $q=4$, a model where the regression coefficients of
the first and second outcomes are equal ($\bm\beta_{1} =
\bm\beta_{2}$), while the coefficients of outcomes three and four are
unrestricted, can be specified as:
% #
% # coef.constraints <- c(1, 1, 2, 3)
% # names(coef.constraints) <- paste0("Y", 1:4)
% # coef.constraints
% # @
\begin{verbatim}
coef.constraints = c(1, 1, 2, 3)
\end{verbatim}
A more flexible framework allows the user to specify such constraints
for each of the regression coefficients of the $p$~covariates, not
only for the whole vector. Such constraints will be specified by means
of a matrix of dimension $q\times p$, where each column specifies
constraints for one of the $p$~covariates in the same way as presented
above. Moreover, a value of \code{NA} indicates that the corresponding
coefficient is fixed (as we will show below) and should not be
estimated.

Let us assume we are dealing with seven covariates and four outcomes
(as in the example presented in Section~\ref{subsect:example1}) and
the following constraints on the regression coefficients:
\begin{itemize}
\item[-] $\beta_{12} = \beta_{22} = \beta_{32}$;
\item[-] $\beta_{13} = 0$, $\beta_{23} = 0$, $\beta_{33} = 0$;
\item[-] $\beta_{14} = \beta_{24} = \beta_{34}$, $\beta_{44} = 0$;
\item[-] $\beta_{15} = \beta_{25} = \beta_{35}=\beta_{45} = 2$.
\end{itemize}

These restrictions on the regression coefficients are imposed by:
\begin{verbatim}
coef.constraints = cbind(c(1, 2, 3, 4), c(1, 1, 1, 2), c(NA, NA, NA, 1),
  c(1, 1, 1, NA), c(NA, NA, NA, NA), c(1, 2, 3, 4), c(1, 2, 3, 4))
\end{verbatim}
Specific values of coefficients can be fixed
through the \code{coef.values} argument, as we will show in the
following.
\subsubsection{Fixing coefficient values}
In addition, specific values on the regression coefficients can be set in
the $q\times p$ matrix \code{coef.values}. Again each column
corresponds to the regression coefficients of one covariate. This
feature is to be used if some of the covariates have known slopes, but
also for excluding covariates from the mean model of some of the
outcomes (by fixing the regression coefficient to zero).

By default, if no \code{coef.values} are passed by the user, all the
regression coefficients which receive an \code{NA} in
\code{coef.constraints} will be set to zero. \code{NA} in the
\code{coef.values} matrix indicates the regression coefficient ought
to be estimated.  Setting \code{coef.values} in accordance with the
\code{coef.constraints} from above:
% # <<
% # coef.values <- cbind(c(NA, NA, NA, NA),
% #                      c(NA, NA, NA, NA),
% #                      c(0, 0, 0, NA),
% #                      c(NA, NA, NA, 0),
% #                      c(2, 2, 2, 2),
% #                      c(NA, NA, NA, NA),
% #                      c(NA, NA, NA, NA))
% # rownames(coef.values) <- paste0("Y", 1:4)
% # colnames(coef.values) <- paste0("X", 1:7)
% # coef.values
% # @
\begin{verbatim}
coef.values = cbind(c(NA, NA, NA, NA), c(NA, NA, NA, NA), c(0, 0, 0, NA),
  c(NA, NA, NA, 0), c(2, 2, 2, 2), c(NA, NA, NA, NA), c(NA, NA, NA, NA))
\end{verbatim}
gives the following equations for the latent scores:
\begin{alignat*}{8}
\widetilde Y_{i1} &= \beta_{11} x_{i1} &&+ \beta_{12} x_{i2} &&+ \beta_{14} x_{i4} &+ 2 x_{i5} &&+ \beta_{16} x_{i6} &&+ \beta_{17} x_{i7} ,\\
\widetilde Y_{i2} &= \beta_{21} x_{i1} &&+ \beta_{12} x_{i2} &&+  \beta_{14} x_{i4} &+ 2 x_{i5} &&+ \beta_{26} x_{i6} &&+ \beta_{27} x_{i7} ,\\
\widetilde Y_{i3} &= \beta_{31} x_{i1} &&+ \beta_{12} x_{i2} &&+ \beta_{14}  x_{i4} &+ 2 x_{i5} &&+ \beta_{36} x_{i6} &&+ \beta_{37} x_{i7},\\
\widetilde Y_{i4} &= \beta_{41} x_{i1} &&+ \beta_{42} x_{i2} &+ \beta_{43} x_{i3} &&+2 x_{i5} &&+ \beta_{46} x_{i6} &&+ \beta_{47} x_{i7}.
\end{alignat*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Constraints on category-specific coefficients}
If the parallel regression or proportional odds assumption ought to be relaxed,
the constraint design of package \pkg{VGAM} can be employed.
Let us consider the model from above and further assume outcome one, two and three are ordinal variables with four categories and outcome four has three categories. For illustration purposes we now relax the parallel regression assumption partially for covariates \code{X1} and \code{X3} in the following way:
\begin{itemize}
\item $\beta_{11, 1} = \beta_{11, 2} \neq \beta_{11,  3}$;

\item $\beta_{43, 1} \neq \beta_{43, 2}$,
\end{itemize}
where $\beta_{jk,  r}$ denotes the regression coefficient of covariate~$k$ in the linear predictor of the $r$-th cumulative probit or logit for measurement~$j$. By the first restriction,  for the first outcome two regression coefficients are employed for covariate \code{X1}: one coefficient $\beta_{11,  1} = \beta_{11, 2}$ for the first two linear predictors and one coefficient $\beta_{11,  3}$ for the third linear predictor. Covariate \code{X3} only appears in the model for the fourth  outcome. For each of the two linear predictors a different regression coefficient is estimated: $\beta_{43, 1}$ and $\beta_{43, 2}$. Note also that some of the coefficients are constant across outcomes. Moreover the coefficient of covariate \code{X5} is fixed to two by setting an offset in the model \code{formula}:
\begin{verbatim}
formula = 0 + X1 + X2 + X3 + X4 + offset(2 * X5) + X6 + X7
\end{verbatim}

The constraints are set up as a named list where the names correspond to the names of all covariates in the \code{model.matrix}. The number of rows is equal to the total number of linear predictors $\sum_j (K_j - 1)$ of the ordered responses, in the example above $3+3+3+2=11$ rows. The number of columns represents the number of parameters to be estimated for each covariate.

\begin{verbatim}
coef.constraints = list(X1 = cbind(c(1, 1, rep(0, 9)), c(0, 0, 1, rep(0, 8)),
  c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 9), rep(1, 2))),
  X2 = cbind(c(rep(1, 9), rep(0, 2)), c(rep(0, 9), rep(0, 2))),
  X3 = cbind(c(rep(0, 9), 1, 0), c(rep(0, 9), 0, 1)),
  X4 = matrix(c(rep(1, 9), rep(0, 2))), X6 = cbind(c(rep(1, 3), rep(0, 8)),
    c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 6), rep(1, 3), rep(0, 2)),
    c(rep(0, 9), rep(1, 2))), X7 = cbind(c(rep(1, 3), rep(0, 8)),
    c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 6), rep(1, 3), rep(0, 2)),
    c(rep(0, 9), rep(1, 2))))
\end{verbatim}
For more details we refer the reader to the documentation of the \pkg{VGAM} package.

%
%
%
%
%
%
%
% If the parallel regression or proportional odds assumption ought to be relaxed,
% the constraint design of package \pkg{VGAM} can be employed.
% Let us consider the model from above and further assume outcome one, two and three are ordinal variables with 3 categories and outcome four has two categories. For illustration purposes we now relax the parallel regression assumption partially for covariates \code{X1} and \code{X3} in the following way:
% \begin{itemize}
% \item $\beta_{11,\text{cat}_1} = \beta_{11,\text{cat}_2} \neq \beta_{11, \text{cat}_3}$;
%
% \item $\beta_{43,\text{cat}_1} \neq \beta_{43,\text{cat}_2}$,
% \end{itemize}
% where $\beta_{jk, \text{cat}_r}$ denotes the regression coefficient for multiple measurement~$j$, covariate~$k$ and category~$r$. By the first restriction,  for the first outcome two regression coefficients are employed for covariate \code{X1}: one coefficient $\beta_{11, \text{cat}_1} = \beta_{11,\text{cat}_2}$ for the first two categories and one coefficient $\beta_{11, \text{cat}_3}$ for the third category. Covariate \code{X3} only appears the model for the fourth  outcome. For each of the two categories a different regression coefficient is estimated: $\beta_{43,\text{cat}_1}$ and $\beta_{43,\text{cat}_2}$. Note also that some of the coefficients are constant across outcomes. Moreover the coefficient of covariate \code{X5} is fixed to two by setting an offset in the model \code{formula}:
% \begin{verbatim}
% formula = 0 + X1 + X2 + X3 + X4 + offset(2 * X5) + X6 + X7
% \end{verbatim}
%
% The constraints are set up as a named list where the names correspond to the names of all covariates in the \code{model.matrix}. The number of rows is equal to the total number of classes of the ordered responses, in the example above $3+3+3+2=11$ rows. The number of columns represents the number of parameters to be estimated for each covariate.
%
% \begin{verbatim}
% coef.constraints = list(X1 = cbind(c(1, 1, rep(0, 9)), c(0, 0, 1, rep(0, 8)),
%   c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 9), rep(1, 2))),
%   X2 = cbind(c(rep(1, 9), rep(0, 2)), c(rep(0, 9), rep(0, 2))),
%   X3 = cbind(c(rep(0, 9), 1, 0), c(rep(0, 9), 0, 1)),
%   X4 = matrix(c(rep(1, 9), rep(0, 2))), X6 = cbind(c(rep(1, 3), rep(0, 8)),
%     c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 6), rep(1, 3), rep(0, 2)),
%     c(rep(0, 9), rep(1, 2))), X7 = cbind(c(rep(1, 3), rep(0, 8)),
%     c(rep(0, 3), rep(1, 3), rep(0, 5)), c(rep(0, 6), rep(1, 3), rep(0, 2)),
%     c(rep(0, 9), rep(1, 2))))
% \end{verbatim}
% For more details we refer the reader to the documentation of the \pkg{VGAM} package.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Interaction terms and categorical covariates}
When constraints on the regression coefficients should be specified in
models with interaction terms or categorical covariates, the
\code{coef.constraints} matrix has to be constructed appropriately. If
the order of the terms in the covariate matrix is not clear to the
user, it is helpful to call the function \code{names_constraints()}
before constructing the \code{coef.constraints} and \code{coef.values}
matrices. Consider the following \code{data}:
<<echo = FALSE>>=
n<-3
data <- cbind.data.frame("Y" = ordered(c("A", "B", "A")),
  X1 = rnorm(n), X2 = rnorm(n),
  X3 = as.factor(c("c1", "c2", "c3")),
  X4 = rnorm(n), X5 = rnorm(n))
data
@
where \code{X3} is a categorical variable with three categories \code{c1}, \code{c2} and \code{c3}.
The names of each column in the covariate matrix can be obtained by:
<<>>=
formula <- Y ~ 1 + X1 : X2 + X3 + X4 + X3 * X5
names_constraints(formula, data)
@
or, without intercept:
<<>>=
formula <- Y ~ 0 + X1 : X2 + X3 + X4 + X3 * X5
names_constraints(formula, data)
@
This should be used when setting up the coefficient constraints. Please note that category \code{c1} for \code{X3} is taken as baseline for both cases.
%%
%
% \begin{alignat*}{7}
% \widetilde Y^{(\text{cat}_1,\text{cat}_2)}_{i1} &= \beta^{(\text{cat}_1,\text{cat}_2)}_{11} x_{i1} &&+ \beta_{12} x_{i2} &&+ \beta_{14} x_{i4} &+ 2 x_{i5} &&+ \beta_{16} x_{i6}  ,\\
% \widetilde Y^{(\text{cat}_3)}_{i1} &= \beta^{(\text{cat}_3)}_{11} x_{i1} &&+ \beta_{12} x_{i2} &&+ \beta_{14} x_{i4} &+ 2 x_{i5} &&+ \beta_{16} x_{i6},\\
% \widetilde Y_{i2} &= \beta_{21} x_{i1} &&+ \beta_{12} x_{i2} &&+  \beta_{14} x_{i4} &+ 2 x_{i5} &&+ \beta_{26} x_{i6} ,\\
% \widetilde Y^{(\text{cat}_1)}_{i3} &= \beta_{31} x_{i1} &&+ \beta_{32} x_{i2} &+ \beta^{(\text{cat}_1)}_{33} x_{i3} &&+2 x_{i5} &&+ \beta_{36} x_{i6} \\
% \widetilde Y^{(\text{cat}_2)}_{i3} &= \beta_{31} x_{i1} &&+ \beta_{32} x_{i2} &+ \beta^{(\text{cat}_2)}_{33} x_{i3} &&+2 x_{i5} &&+ \beta_{36} x_{i6}.
% \end{alignat*}
% This is, for the first outcome, two regression coefficients are employed for covariate \code{X1}: one coefficient $\beta^{(\text{cat}_1,\text{cat}_2)}_{11}$ for the first two categories and one coefficient $\beta^{(\text{cat}_3)}_{11}$ for the third category. Covariate \code{X3} only appears the model for the third  outcome. For each of the two categories a different regression coefficient is estimated: $\beta^{(\text{cat}_1)}_{33}$ and $\beta^{(\text{cat}_2)}_{33}$. Note also that some of the coefficients are constant across outcomes. Moreover the coefficient of covariate \code{X5} is fixed to two.
%
% The constraints are set up as a named list where the names correspond to the names of the covariates. The number of rows is equal to the total number of classes of the ordered responses, in the example above $3+3+2=8$ rows. The number of columns represents the number of parameters to be estimated for each covariate.
% <<eval = F>>=
% coef.constraints = list(
%   X1 = cbind(c(1, 1, rep(0, 6)), c(0, 0, 1, rep(0, 5)), c(0, 0, 0, 1, 1, 1, 0, 0),
%     c(rep(0, 6),1,1)), X2 = cbind(c(rep(1, 6), 0, 0), c(rep(0, 6), 1, 1)),
%   X3 = cbind(c(rep(0, 6), 1, 0), c(rep(0, 6), 0, 1)),
%   X4 = matrix(c(rep(1, 6), rep(0, 2))),
%   X6 = cbind(c(1,1,1,rep(0, 5)), c(0,0,0,1,1,1,0,0), c(rep(0,6), 1,1)))
% @
% The coefficient for covariate \code{X5} fixed to two can be specified in the formula as:
% \begin{verbatim}
% formula = 0 + X1 + X2 + X3 + X4 + offset(2 * X5) + X6
% \end{verbatim}
%
% For more details we refer the reader to the documentation of the \pkg{VGAM} package.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\subsection{Additional arguments}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection[weights]{\code{weights}}
\sloppy
Weights on each subject~$i$ can be chosen in a way that they are
constant across multiple measurements. Weights should be stored in a column of
\code{data}.  The column name of the weights in \code{data} should be
passed to this argument \code{weights}. Negative weights are not allowed.
%
\subsubsection[offset]{\code{offset}}
If offsets are not specified in the model \code{formula}, a list with a vector of offsets for each multiple measurement can be passed.
%
\subsubsection[scale]{\code{scale}}
If \code{scale = TRUE}, the continuous covariates are standardized by substracting the mean and dividing by the standard deviation. This operation is performed for each repeated measurement before fitting.
%
\subsubsection[solver]{\code{solver}}
All general purpose optimizers of the \proglang{R} package
\pkg{optimx} \citep{optimx1, optimx2} can be used for maximization of the composite
log-likelihood. These are \code{"Nelder-Mead", "BFGS", "CG",
  "L-BFGS-B", "nlm", "nlminb", "spg", "ucminf", "newuoa", "bobyqa",
  "nmkb", "hjkb", "Rcgmin"} and \code{"Rvmmin"}. The default is the \code{"newuoa"} solver.
  The \code{"BFGS"} solver performs well in terms of computation time, but it suffers from convergence problems, especially for the \code{mvlogit()} link.
  Moreover, if the user desires a specific solver which is not implemented in the \proglang{R} package \pkg{optimx}, other applicable solvers can be used by using a wrapper function with arguments \code{starting.values}, \code{objFun}, \code{control} of the following form:
\begin{verbatim}
solver = function(starting.values, objFun, control) {
  optRes <- solver.function(...)
  list(optpar = optRes$optpar, objvalue = optRes$objvalue)
}
\end{verbatim}
The \code{solver.function()} should return a list of two elements \code{optpar} and \code{objvalue}.
The element \code{optpar} should be a vector of length equal to number of parameters to optimize containing the estimated parameters, while the element \code{objvalue} should contain the value of the objective function after the optimization procedure.
%
\subsubsection[Standard errors]{\code{se}}
If \code{se = TRUE} standard errors are computed using the Godambe
information matrix (see Section~\ref{sect:estimation}).
%
\subsubsection[Starting values]{\code{start.values}}
A list of starting values for threshold as well as regression
coefficients can be passed by the argument \code{start.values}. This
list contains a list (with a vector of starting values for each
dimension) \code{theta} of all flexible threshold parameters and a
list \code{beta} of all flexible regression parameters. All fixed
values need to be excluded and in case of constraints on a whole
dimension (e.g., \code{threshold.constraints = c(1,1,2,3)} or
\code{coef.constraints = c(1,1,2,3)}), the element can be either
skipped or a vector of length zero can be set. Starting values for
Example 1 in Section~\ref{sect:examples} are for example:
\begin{verbatim}
start.values = list(theta = list(c(-3, -1, 0, 0.5, 2.5),
  c(-3, -1, 0, 0.5, 2, 3.5), c(0)), beta = list(c(0.05, -0.05, -0.8, 1, 0.2),
    c(-0.5, 0.2), c(-0.3, 0.3), c(0.5, -1.1, 0.7, 0.3, -1.2)))
\end{verbatim}
Starting values for the error structure parameters cannot be specified
by the user because of the complexity of the different error structure
parameterizations and transformations.
\subsubsection[contrasts]{\code{contrasts}}
\code{contrasts} can be specified by a named list as in the argument \code{contrasts.arg} in \code{model.matrix.default()}.
\subsubsection[PL.lag]{\code{PL.lag}}
In longitudinal studies, where $q_i$ is possibly large, the pairwise
likelihood estimation can be time consuming as it is build from all
two dimensional combinations of $j\in J_i$. To overcome this
difficulty, one can construct the likelihood using only the bivariate
probabilities for pairs of observations less than $lag$ in ``time
units'' apart. A similar approach was proposed by
\cite{Varin09}. Assuming that, for each subject~$i$, we have a
time-series of consecutive ordinal observations, the $i$-th component
of the pairwise likelihood has the following form:
\begin{align*}
p\ell^{lag}_i(\bm\delta)=
 w_i  \left[ \sum_{k = 1}^{q_i-1} \sum_{l =
k+1}^{q_i}\mathbbm{1}_{\{|t_{k}-t_{l}|\leq lag\}}  \log \Prob(Y_{ik} = r_{ik},
Y_{il} = r_{il}) \right].
\end{align*}
The $lag$ can be fixed by the argument \code{PL.lag}
and it can only be used along with \code{error.structure = cor_ar1()}. The use of this argument is however not recommended if there are missing observations in the time series, i.e., if the ordinal variables are not observed in consecutive years. Moreover, one should also proceed with care if the observations are not missing at random.
%
\subsubsection[Control]{\code{control}}
A list of control arguments that are passed to the function \code{optimx()} or to the user-specific \code{solver.function()}. For further details see \cite{optimx1}.
%
\subsection[clearname]{Methods for class \code{`mvord'}}
Several methods are implemented for the class \code{`mvord'}. These
methods include a \code{summary()} and a \code{print()} function to
display the estimation results, a \code{coef()} function to extract
the regression coefficients, a \code{thresholds()} function to extract
the threshold coefficients and a function \code{get_error_struct()} to
extract the estimated parameters of the correlation/covariance
structure of the errors. The pairwise log-likelihood can be extracted
by the function \code{logLik()}, function \code{vcov()} extracts the
variance-covariance matrix of the parameters and \code{nobs()}
provides the number of subjects. In addition, a \code{predict()}
function, a \code{marginal.predict()} function and a \code{get.prob()}
function are available to predict probabilities, cumulative
probabilities, categories and linear predictors for the fitted model.
Other standard methods such as \code{terms()}, \code{model.matrix()}
and \code{fitted()} are also available. Functions \code{AIC()} and
\code{BIC()} can be used to extract the information criteria CLAIC and
CLBIC.
%
\subsection{Output}
The function \code{mvord()} returns an object of class
\code{`mvord'},
which is a list containing the following
components:

\begin{tabularx}{\textwidth}{lX}
   \code{beta} & a named \code{matrix} of regression coefficients\\
   \code{theta} & a named \code{list} of threshold parameters \\
\code{error.struct} & an object of class \code{`error\_struct'}\\
   \code{sebeta} & a named \code{matrix} of standard errors of the regression coefficients\\
   \code{setheta} & a named \code{list} of standard errors of the threshold parameters\\
\code{seerror.struct} & a \code{vector} of standard errors for the parameters of the error structure\\
   \code{rho} & a \code{list} of all objects that are used in \code{mvord()}
   \end{tabularx}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection[clearname]{Implementation \code{mvord2()}}
Additionally, a second function \code{mvord2()} is implemented, for the
setting where the covariates do not vary between the multiple
measurements ($\bm x_{i1}=\ldots= \bm x_{iq}$):
\begin{verbatim}
mvord2(formula, data, error.structure = cor_general(~ 1), link = mvprobit()
  coef.constraints = NULL, coef.values = NULL, threshold.constraints = NULL,
  threshold.values = NULL, weights = NULL, se = TRUE, offset = NULL,
  scale = FALSE, start.values = NULL, solver = "newuoa", PL.lag = NULL,
  control = list(maxit = 200000, trace = 1, kkt = FALSE))
\end{verbatim}
This function uses a slightly simplified data structure, where the multiple ordinal observations as well as the covariates
are stored as columns in a \code{data.frame}. Each subject~$i$ corresponds to one row of the data
frame, where all outcomes $Y_{i1}, \dots, Y_{iq}$ (with missing observations set to \code{NA})
and all the covariates $x_{i1}, \dots, x_{ip}$ are stored in different columns.
Each outcome must be an ordered factor.
In order to specify the model we use a
multivariate \code{formula} object of the form:
\begin{verbatim}
formula = cbind(Y1, ..., Yq) ~ 0 + X1 + ... + Xp
\end{verbatim}
The \code{error.structure} and the constraints on the regression and
threshold parameters are set in analogy to \code{mvord()}. However,
the ordering of the responses is given by the ordering in the model
\code{formula}. In addition, the \code{link}, subject \code{weights},
\code{se} and the \code{solver} are chosen in the same way as in
\code{mvord()}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DATA DESCRIPTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Description of the illustrative data sets}\label{sect:data}
As previously mentioned, the motivation of the \pkg{mvord} package
lies in a credit risk application, where ordinal credit ratings from
different rating sources are available for a sample of firms over
several years. In this application, yearly credit ratings from the big
three credit rating agencies Standard \& Poor's, Moody's and Fitch
were collected for a panel of publicly traded US companies. Moreover,
accounting, stock price and default information on a yearly basis
have been obtained for each firm-year in the sample.

Before proceeding with the multivariate analysis of these credit
ratings, the following aspects of the problem at hand are worth
mentioning. Firstly, a feature of the sample is the rather high
number of missing values. Not all firms in the sample are rated by all
three raters at the same point in time. In line with the market share
the raters hold in the US ratings market, the number of rated firms
per year differs between the three raters.  Secondly, the raters can
employ different rating scales in their rating methodology, e.g.,
different labeling of their ordinal rating classes or a different
number of rating classes. A perhaps more fundamental difference
between the rating scales appears when the raters employ different
rating ``philosophies'', i.e., when the central credit risk measure
employed in assessing creditworthiness differs among raters (e.g.,
probability of default versus recovery in case of default).
This also has implications for the covariates used in the rating
models, as the inclusion or importance of various covariates can
differ among the raters. Even though the rating agencies do not fully
disclose their rating process, most likely for keeping their
competitive advantage, they do publish reports on a regular basis
where they state what are the main variables considered in the rating
process.

These data have been analyzed in \cite{epubwu5389}, where a
multivariate model of corporate credit ratings has been proposed.  Unfortunately, this original
data set is not freely redistributable. Therefore, we resort to the
simulation of illustrative data sets by taking into consideration key
features of the original data such as the ones mentioned above. The simulated data sets contain ratings from four different rating
sources and firm-specific covariates for a panel of $\Sexpr{length(unique(data_cr_mvord2$firm_id))}$ firms over
a sample period of ten years.

%% Covariates
\subsection{Covariates}
The covariates employed correspond to firm-level and market financial ratios which measure different aspects of a firm's financial health:

\begin{itemize}
\item \code{ICR} -- interest coverage ratio, which measures how well the interest expenses can be covered from the free operating cash-flow of a company; % R5)
\item \code{LR} -- liquidity ratio relating the cash held by a company to the current liabilities; % R8 -- Cash to current liabilities
\item \code{LEV1} -- leverage ratio relating debt to earnings before interest and taxes; %R24 Debt to EBITDA
\item \code{LEV2} -- leverage ratio measuring the percentage of debt in the long-term capital of a firm; % R31
\item \code{PR} -- profitability ratio measuring return on capital;% R39
%\item \code{CFR} -- cash-flow to debt ratio; % R44 -
\item \code{lRSIZE} -- relative size of the company in the market (in logarithm);  %RSIZE
\item \code{lSYSR} -- a measure of systematic risk (in logarithm);%SIGMA
\item \code{BSEC} -- the business sector of a firm.
\end{itemize}
More information on these covariates can be found in e.g., \cite{campbell2008search} or \cite{sp}. The business sector variable \code{BSEC} has been kept as in the
original data set.  For the continuous covariates we resort to
simulation. As the distribution of the continuous covariates in
the original data set is rather homogeneous across the ten years, we fit a
distribution to each covariate over the whole sample period using the
function \code{fitdistr()} of the \pkg{MASS} package. The best fitting
distribution has been chosen by AIC.
%% In line with previous results, most of the ratios have
%% heavy tails.
The parameters and distributions used for the simulation of the
covariates are given in Table~\ref{tab:covariates}.
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Covariate & Bound & Distribution & Parameters \\
  \hline
\code{ICR} & unbounded & $t(m, s, df)$ & $m = 2.301$, $s = 2.120$, $df = 1.163$\\
  \code{LR} & positive & log-$t(m, s, df)$ & $m = -1.472$, $s = 0.992$, $df = 2.753$ \\
  \code{LEV1}  & positive & log-$t(m, s, df)$& $m = 1.008$, $s = 0.596$, $df = 1.245$ \\

  \code{LEV2} & positive & log-$t(m, s, df)$& $m = -0.777$, $s = 0.443$, $df = 1.289$ \\

  \code{PR}  & unbounded & $t(m, s, df)$ & $m = 0.12$, $s = 0.08$, $df = 4$ \\
  \code{lRSIZE} & unbounded &$\mathcal{N}(m,s)$ & $m = -9.131$, $s = 1.475$ \\

  \code{lSYSR} & unbounded & $\mathcal{N}(m,s)$ & $m = -3.877$, $s = 0.439$ \\
   \hline
\end{tabular}
\caption{Distributions used for the simulation of the covariates.}
\label{tab:covariates}
\end{table}

After simulating the covariates from the distributions presented in Table~\ref{tab:covariates}, we winsorize the simulated ratios at 10\% for unbounded ratios and 5\% for bounded ratios (by setting for all ratios the values above the 95\%-quantile to the 95\%-quantile; for the unbounded ratios we also winsorize from below, by setting for the values below the 5\%-quantile to the 5\%-quantile). This is standard practice in the literature on credit risk estimation based on financial ratios for reducing the influence of outliers \citep[e.g.,][]{campbell2008search}.
%% Cross-section
\subsection{Cross-sectional data set}\label{subsect:crossdata}
For the $\Sexpr{length(unique(data_cr_mvord2[,'firm_id']))}$ firms in
the first year of the sample, we generate ratings from four rating
sources.  We name the four rating sources \code{R1}, \code{R2},
\code{R3} and \code{R4}. Raters \code{R1}, \code{R2} and \code{R3}
correspond to the three rating agencies in the original
application. The fourth rating source \code{R4} corresponds to an
indicator variable constructed from the original data set, which
indicates if the company will be in an investment or speculative grade
within one year from the rating observations. In the original data
set, this indicator is computed by checking the one-year ahead rating
for each of the three raters. We then follow the US regulation stating
that, if two ratings for a company are available, one investment grade
and one speculative grade, the company is categorized as speculative
grade, whereas if three ratings are available the majority rule
applies \citep[][]{JOFI:JOFI1709}.

Similarly to the original data, we assume
that the rating methodology of the raters is heterogeneous. Raters
\code{R1} and \code{R2} employ a similar credit risk measure when
assigning firms on a six-point ordinal scale (which we label, from
best to worst, \code{A}, \code{B}, \code{C}, \code{D}, \code{E},
\code{F}), rater \code{R3} uses a seven-point scale and different
labeling (from best to worst, \code{G}, \code{H}, \code{I}, \code{J},
\code{K}, \code{L}, \code{M}) while the \code{R4} ratings differentiate on a two
point scale between speculative grade and investment grade firms (from best to worst, \code{N}, \code{O}).

We simulate the latent scores $\widetilde Y_{ij}$ from the model in
Equation~\ref{eqn:model} with multivariate logit link and a constant general correlation structure.   The
correlation matrix used in the simulation is
\[\Sigma= \left( {\begin{array}{cccc}
   1 & 0.9 & 0.7 & 0.5 \\
   0.9 & 1 & 0.7 &0.6\\
   0.7 & 0.7 & 1 & 0.8 \\
   0.5 & 0.6 &0.8 & 1\\
  \end{array} } \right).
\]
We use the following coefficients:
\begin{itemize}
\item $\beta_{0j}=0$ $\forall j\in\{\text{\code{R1}, \code{R2}, \code{R3}, \code{R4}}\}$;
\item $\bm \beta_{\text{\code{R1}}} = (0.3, 0, -0.2,    -1.3,  0.5, 0.2, 0)^\top$;
\item $\bm \beta_{\text{\code{R2}}} = (0,   0, -0.2,    -1,    0.5, 0.3, 0)^\top$;
\item $\bm \beta_{\text{\code{R3}}} = (0.3, 0, -0.2,    -1.4,  0.5,   0.4, 0)^\top$;
\item $\bm \beta_{\text{\code{R4}}} = (0,   -0.3,  0,  -1.6,  1.9, 0.1, -0.2)^\top$.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The threshold parameters are chosen such that the simulated
distribution of ratings (the proportion of observations falling into
each rating class) is similar to the distribution of the ratings in
the original data set.

As raters \code{R1} and \code{R2} employ similar rating
methodologies and use the same labeling of their rating classes, we
set $\bm\theta_{\text{\code{R1}}}=\bm\theta_{\text{\code{R2}}}=(-4.5,
-3, -1, 1, 3.5)^\top$. The thresholds for rater~\code{R3} are:
$\bm\theta_{\text{\code{R3}}}=(-6,-3 ,0 ,1, 2, 4)^\top$ and the
threshold for rater~\code{R4} is $\theta_{\text{\code{R4}}}=-1$.

The panel of ratings in the original data set is highly unbalanced, as
not all firms receive ratings from all four sources. We therefore keep
the missingness pattern and remove the simulated ratings that
correspond to missing observations in the original data.  For
rater~\code{R1} we remove
$\Sexpr{round(sum(is.na(data_cr_mvord2[,'R1']))/NROW(data_cr_mvord2)*100,2)}$\%,
for rater~\code{R2}
$\Sexpr{round(sum(is.na(data_cr_mvord2[,'R2']))/NROW(data_cr_mvord2)*100,2)}$\%,
and for rater~\code{R3}
$\Sexpr{round(sum(is.na(data_cr_mvord2[,'R3']))/NROW(data_cr_mvord2)*100,2)}$\%
of the observations.

The cross-sectional data set is available in the \pkg{mvord} package as \code{data_cr_mvord} (in long format) or  \code{data_cr_mvord2} (in short format). More details on the structure of this data set will be provided in Section~\ref{sect:examples}.
%% Longitudinal
\subsection{Longitudinal data set}
For the whole panel of firms we simulate ratings
corresponding to rater~\code{R1} over a period of ten years (\code{year1}, \code{year2}, $\dots$, \code{year10}). From the initial
$\Sexpr{NROW(data_cr_mvord2)}$ firms which are in the sample in the
first year,
$\Sexpr{length(unique((data_cr_panel[,'firm_id'])))}$ are rated by
rater~\code{R1}.

We simulate the latent scores from the model in Equation~\ref{eqn:model} with probit link and business sector-specific $AR(1)$ correlation structure. The dependence of the correlation structure on the business sector is motivated by the fact that in some sectors such as manufacturing ratings tend to be more ``sticky'', i.e., do not change often over the years, while in more volatile sectors like IT there is less ``stickiness'' in the ratings. The sector specific correlations are: $\rho_{\text{\code{BSEC1}}}=0.9$,
 $\rho_{\text{\code{BSEC2}}}=0.7$,
 $\rho_{\text{\code{BSEC3}}}=0.9$,
 $\rho_{\text{\code{BSEC4}}}=0.9$,
 $\rho_{\text{\code{BSEC5}}}=0.9$,
 $\rho_{\text{\code{BSEC6}}}=0.7$,
 $\rho_{\text{\code{BSEC7}}}=0.5$,
 $\rho_{\text{\code{BSEC8}}}=0.6$.

We use no intercepts and two sets of regression coefficients:
\begin{itemize}
\item $\bm\beta_{\text{\code{year1}}}=\dots= \bm\beta_{\text{\code{year5}}}=(0.17,0,-0.11,-0.72,0.28,0.11,0)^\top$;
\item $\bm\beta_{\text{\code{year6}}}=\dots= \bm\beta_{\text{\code{year10}}}=(0.08,-0.2,-0.3,-1.6,0.4,0.5,-0.2)^\top$.
\end{itemize}
Such a ``break'' in the coefficients could be explained by e.g., a change of regimes. Moreover, we assume the rating methodology of rater~\code{R1} does not change
over time, hence we use one set of threshold parameters for all years: $\bm\theta_{\text{\code{year1}}}=\dots=\bm\theta_{\text{\code{year10}}}=( -4, -1, 0, 1, 3)^\top$. The longitudinal panel is unbalanced, as firms can leave the sample
for various reasons such as default or mergers and acquisitions. In addition, credit ratings can be withdrawn. We remove the simulated ratings that
correspond to missing firm-year observations in the original data. The longitudinal panel of \code{R1} ratings together with the covariates is available in \code{data_cr_panel}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Examples}\label{sect:examples}
We use the credit risk data sets presented in Section~\ref{sect:data}
to illustrate some of the features of package \pkg{mvord}.
%
\subsection{Example 1 -- Ratings assigned by multiple raters to a cross-section of firms}\label{subsect:example1}
The first example presents a multivariate ordinal regression model with logit link and a general correlation error structure (\code{cor_general(~ 1)}).
The simulated data set contains the credit risk measure \code{rating} (ratings assigned by raters \code{R1}, \code{R2}, \code{R3} and \code{R4}) and 8 covariates for a cross-section of \Sexpr{length(unique(data_cr_mvord[,"firm_id"]))} firms.
The number of firm-ratings is \Sexpr{NROW(data_cr_mvord)}.
<<>>=
head(data_cr_mvord, n = 3)
str(data_cr_mvord, vec.len = 3)
@
The distribution of the ratings classes for the four raters is:
<<>>=
by(data_cr_mvord,  data_cr_mvord$rater_id,
    function(x) table(x$rating))
@
We include 7 financial ratios as covariates in a model without intercept by the \code{formula}:
\begin{verbatim}
formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR
\end{verbatim}
The subject index $i$ is stored in the column \code{firm_id} and the multiple measurement index $j$, which indicates the rater, is given by \code{rater_id}:
\begin{verbatim}
index = c("firm_id", "rater_id")
\end{verbatim}
An optional vector \code{response.names} is used to specify all the raters to be included in the model. The ordering of this vector is essential when constraints on the parameter set want to be imposed:
\begin{verbatim}
response.names = c("R1", "R2", "R3", "R4")
\end{verbatim}
Due to the fact that the categories differ across raters we specify the \code{response.levels} by:
\begin{verbatim}
response.levels = list(rev(LETTERS[1:6]), rev(LETTERS[1:6]),
  rev(LETTERS[7:13]), rev(LETTERS[14:15]))
\end{verbatim}
If no \code{response.levels} are passed, the natural ordering is used and could lead to an incorrect labeling.
The rating classes assigned by the raters are here in order from worst to best indicating that lower values of the latent variables indicate lower creditworthiness or increased credit risk.

We fit a model to these data with the following features:
\begin{itemize}
 \item We assume that \code{R1} and \code{R2} use the same rating scale by setting the following constraints on the threshold parameters:
\begin{verbatim}
 threshold.constraints = c(1, 1, 2, 3)
\end{verbatim}
\item We assume that some covariates are equal for some raters. For example, we assume that the coefficient of \code{ICR} is equal for \code{R1} and \code{R3}, or that the coefficients of \code{LEV1} and \code{PR} are the same for the raters \code{R1}, \code{R2} and \code{R3}. In addition, some of the regression coefficients are set to zero like \code{ICR} for \code{R1} and \code{R3}, or \code{lSYSR} for the raters \code{R1}, \code{R2} and \code{R3}. All the constraints above and some additional constraints are performed by the following restrictions on the regression coefficients by using the more flexible method:
\begin{verbatim}
coef.constraints = cbind(c(1, NA, 1, NA), c(NA, NA, NA, 1),
  c(1, 1, 1, NA), c(1, 2, 3, 4), c(1, 1, 1, 4), c(1, 2, 3, 4),
  c(NA, NA, NA, 1))
\end{verbatim}
The \code{NA}s in \code{coef.constraints} are fixed to specific values. If no matrix \code{coef.values} is provided, the coefficients are set by default to zero automatically. This automatically generated \code{coef.values} matrix, looks like:
\begin{verbatim}
coef.values = cbind(c(NA, 0, NA, 0), c(0, 0, 0, NA), c(NA, NA, NA, 0),
  c(NA, NA, NA, NA), c(NA, NA, NA, NA), c(NA, NA, NA, NA),
  c(0, 0, 0, NA))
\end{verbatim}
The specified \code{coef.constraints} together with \code{coef.values} give the following model:
\begin{alignat*}{8}
\widetilde Y_{1} &= \beta_{11} \text{\code{ICR}}+ && & \beta_{13}\text{\code{LEV1}}& + &\beta_{14} \text{\code{LEV2}} &+ \beta_{15}\text{\code{PR}}  &+ \beta_{16} \text{\code{lRSIZE}} ,\\
\widetilde Y_{2} &=   && &\beta_{13} \text{\code{LEV1}}&+& \beta_{24} \text{\code{LEV2}}&+ \beta_{15} \text{\code{PR}} &+ \beta_{26} \text{\code{lRSIZE}}  ,\\
 \widetilde Y_{3} &= \beta_{11}\text{\code{ICR}}+  &&& \beta_{13} \text{\code{LEV1}} &+ &\beta_{34}  \text{\code{LEV2}}&+ \beta_{15} \text{\code{PR}} &+ \beta_{36}\text{\code{lRSIZE}} ,\\
 \widetilde Y_{4} &= & \beta_{42} \text{\code{LR}}+ && & &\beta_{44}\text{\code{LEV2}} &+ \beta_{45}\text{\code{PR}}  &+ \beta_{46} \text{\code{lRSIZE}} &+ \beta_{47} \text{\code{lSYSR}}.
\end{alignat*}
%\item allows for different correlation error structures depending on the business sector:
\end{itemize}
As a link function we choose the multivariate logit link:
\begin{verbatim}
link = mvlogit()
\end{verbatim}
For simplicity, we use a general correlation structure which is constant for all subjects:
\begin{verbatim}
error.structure = cor_general(~ 1)
\end{verbatim}
In order to avoid numerical instabilities, we standardize our data for each rater:
<<>>=
covar_names <- c("ICR", "LR", "LEV1", "LEV2", "PR", "lRSIZE", "lSYSR")
 data_cr_mvord_scaled <- do.call("rbind.data.frame",
   by(data_cr_mvord, data_cr_mvord$rater_id,
     function(x){x[, covar_names] <- scale(x[, covar_names]); x}))
@
The estimation can now be performed by the function \code{mvord()}:
<<eval = FALSE>>=
res_cor_logit <- mvord(
  formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR,
  data = data_cr_mvord_scaled, error.structure = cor_general(~ 1),
  link = mvlogit(), index = c("firm_id", "rater_id"),
  response.names = c("R1", "R2", "R3", "R4"),
  response.levels = list(rev(LETTERS[1:6]), rev(LETTERS[1:6]),
    rev(LETTERS[7:13]), rev(LETTERS[14:15])),
  coef.constraints = cbind(c(1, NA, 1, NA), c(NA, NA, NA, 1),
    c(1, 1, 1, NA), c(1, 2, 3, 4),  c(1, 1, 1, 4), c(1, 2, 3, 4),
    c(NA, NA, NA, 1)), threshold.constraints = c(1, 1, 2, 3),
  solver = "newuoa")
@
%
<<echo = FALSE, results = hide>>=
FILE <- "res_cor_logit.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
    res_cor_logit <- mvord(
      formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR,
      error.structure = cor_general(~ 1),
      link = mvlogit(),
      data = data_cr_mvord_scaled,
      index = c("firm_id", "rater_id"),
      response.names = c("R1", "R2", "R3", "R4"),
      response.levels = list(rev(LETTERS[1:6]),
                         rev(LETTERS[1:6]),
                         rev(LETTERS[7:13]),
                         rev(LETTERS[14:15])),
      coef.constraints = cbind(c(1, NA, 1, NA),
                               c(NA, NA, NA, 1),
                               c(1, 1, 1, NA),
                               c(1, 2, 3, 4),
                               c(1, 1, 1, 4),
                               c(1, 2, 3, 4),
                               c(NA, NA, NA, 1)),
      threshold.constraints = c(1, 1, 2, 3),
      solver = "newuoa")
  save(res_cor_logit, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }

}
@
%% Then discuss output
The results are displayed by the function \code{summary()}:
<<>>=
summary(res_cor_logit, call = FALSE)
@
If the user fixes thresholds to specific values, the $z$~values and the corresponding $p$~values as well as the significance codes are set to \code{NA}. Another option to display the results is the function \code{print()}:
<<>>=
print(res_cor_logit, call = FALSE)
@
An extended summary, where all thresholds and regression coefficients are shown, even though they are duplicated, can be obtained by:
<<eval = FALSE>>=
summary(res_cor_logit, short = FALSE, call = FALSE)
@
The threshold coefficients can be extracted by the function \code{thresholds()}:
<<>>=
thresholds(res_cor_logit)
@
The regression coefficients are obtained by the function \code{coef()}:
<<>>=
coef(res_cor_logit)
@
The error structure for firm with \code{firm_id} 11 is displayed by the function \code{get_error_struct()}:
<<>>=
get_error_struct(res_cor_logit)[["11"]]
@
%
\subsubsection[clearname]{Fitting the model with the function \code{mvord2()}}
Due to the fact that the covariates do not change across the multiple
measurements (the covariates are firm-specific and do not vary across
raters), we can alternatively fit the model by the function
\code{mvord2()}. In \code{mvord2()}, a slightly different format of \code{data} is used and the ordering of the responses is defined by a
multivariate \code{formula} object. The multiple measurements are stored in
different columns as ordered factors:
<<>>=
head(data_cr_mvord2, n = 3)
str(data_cr_mvord2, vec.len = 2)
@
Again, we standardize the data to avoid numerical instabilities:
<<>>=
data_cr_mvord2[, covar_names] <- scale(data_cr_mvord2[, covar_names])
@
The estimation is performed by calling the function \code{mvord2()}:
<<eval = FALSE>>=
res_cor_logit <- mvord2(
  formula = cbind(R1, R2, R3, R4) ~ 0 + ICR + LR + LEV1 + LEV2 + PR +
    lRSIZE + lSYSR, error.structure = cor_general(~ 1), link = mvlogit(),
  data = data_cr_mvord2, coef.constraints = cbind(c(1, NA, 1, NA),
    c(NA, NA, NA, 1), c(1, 1, 1, NA), c(1, 2, 3, 4), c(1, 1, 1, 4),
    c(1, 2, 3, 4), c(NA, NA, NA, 1)),
  threshold.constraints = c(1, 1, 2, 3))
@
yielding equivalent results to the fit of \code{mvord()}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Example 2 -- Ratings assigned by one rater to a panel of firms}
In a second example we present a longitudinal multivariate ordinal probit regression model with a covariate dependent $AR(1)$ error structure.
The simulated data set contains the credit risk measure \code{rating} (ratings assigned by rater \code{R1}) and 8 covariates for a panel of \Sexpr{length(unique(data_cr_panel[, "firm_id"]))} firms over ten years. The number of firm-year observations is \Sexpr{NROW(data_cr_panel)}:
<<>>=
str(data_cr_panel, vec.len = 3)
head(data_cr_panel, n = 3)
@
The panel is highly unbalanced. The distribution of the number of ratings per firm assigned by rater \code{R1} over the 10 years is given by:
<<>>=
summary(rowSums(with(data_cr_panel, table(firm_id, year))))
@
Per year the number of ratings in the data set decreases:
<<>>=
with(data_cr_panel, table(year))
@
%% first discuss input
We include the 7 financial ratios as covariates in a model without intercept by the \code{formula}:
\begin{verbatim}
formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR
\end{verbatim}
The subject index~$i$ is stored in the column \code{firm_id} while the multiple measurement index~$j$ is given in the column \code{year}:
\begin{verbatim}
index = c("firm_id", "year")
\end{verbatim}
If we wish to estimate the model only for the last eight years of the sample, this can be done by
specifying the names of each ordered response which should enter the model:
\begin{verbatim}
response.names = paste0("year", 3:10)
\end{verbatim}
The rating classes assigned by rater \code{R1} are:
<<>>=
levels(data_cr_panel$rating)
@
with the sixth rating class \code{F} being the worst class and the first rating class \code{A} being the best rating class. We specify the response levels, in the order from worst to best, for each of the eight outcome dimensions through the \code{response.level} argument. Ordering the classes from worst to best indicates that lower values of the latent variables indicate lower creditworthiness or increased credit risk. The rating classes and labels do not change over the eight years:
\begin{verbatim}
response.levels = rep(list(levels(data_cr_panel$rating)),
  length(response.names))
\end{verbatim}
Additionally, the model has the following features:
\begin{itemize}
\item We assume that the rating agencies do not change their methodology over the sample period. This means the threshold parameters are constant over the years. This can be specified through the argument \code{threshold.constraints}:
\begin{verbatim}
threshold.constraints = rep(1, length(response.names))
\end{verbatim}
\item We assume that there is a break-point in the regression coefficients after \code{year5} in the sample. This break-point could correspond to the beginning of a crisis in a real case application.
Hence, we use one set of regression coefficients for years \code{year3}, \code{year4} and \code{year5} and a different set for \code{year6}, \code{year7}, \code{year8}, \code{year9}, \code{year10}.
This can be specified through the argument \code{coef.constraints}:
\begin{verbatim}
coef.constraints = c(rep(1, 3),  rep(2, 5))
\end{verbatim}
\item We allow for different correlation parameters in the $AR(1)$ structure for the different business sectors:
\begin{verbatim}
error.structure = cor_ar1(~ BSEC)
\end{verbatim}
\end{itemize}
As before, we standardize our covariates on a yearly basis:
<<>>=
data_cr_panel_scaled <- do.call("rbind.data.frame",
  by(data_cr_panel, data_cr_panel$year,
    function(x){x[, covar_names] <- scale(x[, covar_names]); x}))
@
The estimation is performed by calling the function \code{mvord()}:
<<eval=F>>=
res_AR1_probit <- mvord(
  formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR,
  error.structure = cor_ar1(~ BSEC), link = mvprobit(),
  index = c("firm_id", "year"), data = data_cr_panel_scaled,
  response.levels = rep(list(levels(data_cr_panel$rating)), 8),
  response.names = paste0("year", 3:10), coef.constraints = c(rep(1, 3),
    rep(2, 5)), threshold.constraints = rep(1, 8), solver = "BFGS")
@
<<echo = FALSE, results = hide>>=
FILE <- "res_AR1_probit.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
    res_AR1_probit <- mvord(
      formula = rating ~ 0 + ICR + LR + LEV1 + LEV2 + PR + lRSIZE + lSYSR,
      error.structure = cor_ar1(~ BSEC),
      link = mvprobit(),
      data = data_cr_panel_scaled,
      index = c("firm_id", "year"),
      response.names = paste0("year", 3:10),
      response.levels = rep(list(levels(data_cr_panel$rating)), 8),
      coef.constraints =  c(rep(1, 3),  rep(2, 5)),
      threshold.constraints = rep(1,8))
  save(res_AR1_probit, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }

}
@
%% Then discuss output
The results are displayed either by the function \code{summary()}:
<<>>=
summary(res_AR1_probit, short = TRUE, call = FALSE, digits = 6)
@
or by the function \code{print()}:
<<>>=
print(res_AR1_probit, call = FALSE, digits = 4)
@
An extended summary, where all thresholds and regression coefficients are shown, even though they are duplicated, can be obtained by:
<<eval = F>>=
summary(res_AR1_probit, short = FALSE, call = FALSE)
@
The threshold coefficients can be extracted by the function \code{thresholds()}:
<<>>=
thresholds(res_AR1_probit)
@
The regression coefficients are obtained by the function \code{coef()}:
<<>>=
coef(res_AR1_probit)
@
The error structure is displayed by the function \code{get_error_struct()}:
<<>>=
get_error_struct(res_AR1_probit)
@
In addition, the correlation parameters $\rho_i$ for each firm are obtained by:
<<>>=
head(get_error_struct(res_AR1_probit, type = "corr"), n = 3)
@
Moreover, the correlation matrices for each specific firm are obtained by:
<<>>=
head(get_error_struct(res_AR1_probit, type = "sigmas"), n = 1)
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Conclusion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion} \label{sect:conclusion}
The present paper is meant to provide a general overview on the
\proglang{R} package \pkg{mvord}, which implements the estimation of
multivariate ordered probit and logit regression models using the
pairwise likelihood approach. Different error structures like a
general correlation and covariance structures, a covariate dependent
equicorrelation structure and an $AR(1)$ structure can be imposed. In
addition, the flexible modeling framework allows imposing constraints
on threshold as well as regression coefficients. Two different data formats
can be used either by applying \code{mvord()} or \code{mvord2()} to estimate the model parameters.

Further research and possible extensions of \pkg{mvord} could be addressed to the
implementation of variable selection procedures in multivariate ordinal regression models
and the inclusion of multivariate semi- or nonparametric ordinal models.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{mvord}
\clearpage
\end{document}
