mvord : An R Package for Fitting Multivariate Ordinal Regression Models

The R package mvord implements composite likelihood estimation in the class of multivariate ordinal regression models with a multivariate probit and a multivariate logit link. A ﬂexible 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 diﬀerent error structures. Heterogeneity in the error structure across the subjects can be accounted for by the package, which allows for covariate dependent error structures. In addition, diﬀerent regression coeﬃcients and threshold parameters for each response are supported. If a reduction of the parameter space is desired, constraints on the threshold as well as on the regression coeﬃcients can be speciﬁed by the user. The proposed multivariate framework is illustrated by means of a credit risk application.


Introduction
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, poor to very good or strongly disagree to 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., modeling 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.).In such a multi-variate setting, models which are able to deal with the correlation in the ordinal outcomes are desired.One possibility is to employ a 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 (see, e.g., Fahrmeir and Tutz 2001).
Several ordinal regression models can be employed for the analysis of ordinal data, with cumulative link models being the most popular ones (e.g., Tutz 2012;Christensen 2019a).Other approaches include continuation-ratio or adjacent-category models (e.g., Agresti 2002Agresti , 2010)).Different packages to analyze and model ordinal data are available in R (R Core Team 2020).For univariate ordinal regression models with fixed effects the function polr() of the MASS package (Venables and Ripley 2002), the function clm() of the ordinal package (Christensen 2019b), which supports scale effects as well as nominal effects, and the function vglm() of the VGAM package (Yee 2010) are available.Another package which accounts for heteroskedasticity is oglmx (Carroll 2018).Package ordinalNet (Wurm, Rathouz, and Hanlon 2017) offers tools for model selection by using an elastic net penalty, whereas package ordinalgmifs (Archer, Hou, Zhou, Ferber, Layne, and Gentry 2014) performs variable selection by using the generalized monotone incremental forward stagewise (GMIFS) method.Moreover, ordinal logistic models can be fitted by the functions lms() and orm() in package rms (Harrell Jr 2019), while ordinal probit models can be fitted by the MCMCoprobit() function in package MCMCpack (Martin, Quinn, and Park 2011) 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 Stata (StataCorp.2018), SAS (SAS Institute Inc. 2018b) or SPSS (IBM Corporation 2017) is provided by Liu (2009).These software packages include the Stata procedure OLOGIT, the SAS procedure PROC LOGISTIC and the SPSS procedure PLUM which perform ordinal logistic regression models.The software procedure PLUM additionally includes other link functions like probit, complementary log-log, cauchit and negative log-log.Ordinal models for multinomial data are available in the SAS package PROC GENMOD, while another implementation of ordinal logistic regression is available in JMP (SAS Institute Inc. 2018a).In Python (Python Software Foundation 2018), package mord (Pedregosa-Izquierdo 2015) implements ordinal regression methods.
While there are sufficient software tools in R which deal with the univariate case, the readyto-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 R packages which support correlated ordinal data.One-dimensional normally distributed random effects in ordinal regression can be handled by the clmm() function of package ordinal (Christensen 2019b).Multiple possibly correlated random effects are implemented in package mixor (Hedeker, Archer, Nordgren, and Gibbons 2018).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 package brms (Bürkner 2017).Multivariate ordinal probit models, where the subject errors are assumed to follow a multivariate normal distribution with a general correlation matrix, can be estimated with package PLordprob (Kenne Pagui and Canale 2018), 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 the package does not account for missing observations in the outcome variable.Polychoric correlations, which are used to measure association among two ordinal outcomes, can be estimated by the polychor() function of package polycor (Fox 2019), where a simple bivariate probit model without covariates is estimated using maximum likelihood estimation.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 nlme (Pinheiro, Bates, and R Core Team 2020), even though models dealing with ordinal data are not supported.
The original motivation for 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 is related to the ability of the firm to meet its financial obligations.In the literature, this latent variable motivation has been used in various credit rating models (e.g., Blume, Lim, and Mackinlay 1998;Afonso, Gomes, and Rother 2009;Alp 2013;Reusens and Croux 2017).This setting is an example of an application where correlated ordinal data arise naturally.On the one hand, multiple ratings assigned by different raters to one firm at the same point in time can be assumed to be correlated.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, there is heterogeneity in the rating methodology.Raters use different labeling as well as a different number of rating classes.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, the data have missing values and are 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 raters at each time point.
The scope of the application of multivariate ordinal regression models reaches far beyond credit risk applications.For example, pain severity studies are a popular setting where repeated ordinal measurements occur.A migraine severity study was employed by Varin and Czado (2010), where patients recorded their pain severity over some time period.In addition to a questionnaire with personal and clinical information, covariates describing the weather conditions were collected.Another application area constitutes the field of customer satisfaction surveys, where questionnaires with ordinal items are often divided into two separate blocks (e.g., Kenne Pagui and Canale 2016).A first block contains questions regarding the general importance of some characteristics of a given service, and a second block relates more to the actual satisfaction on the same characteristics.An analysis of the dependence structure between and within the two blocks is of particular interest.Furthermore, in the presence of multi-rater agreement data, where several raters assign ordinal rankings to different indi-viduals, the influence of covariates on the ratings can be investigated and an analysis and a comparison of the rater behavior can be conducted (e.g., DeYoreo and Kottas 2018).In addition to these few examples mentioned above, the class of multivariate ordinal regression models implemented in mvord (Hirk, Hornik, and Vana 2020b) can be applied to other settings where multiple or repeated ordinal observations occur.This paper discusses package mvord for R which aims at providing a flexible framework for analyzing correlated ordinal data by means of the class of multivariate ordinal regression models and which is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=mvord.In this model class, each of the ordinal responses is modeled as a categorized version of an underlying continuous latent variable which is slotted according to some threshold parameters.On the latent scale we assume a linear model for each of the underlying continuous variables and the existence of a joint distribution for the corresponding error terms.A common choice for this joint distribution is the multivariate normal distribution, which corresponds to the multivariate probit link.We extend the available software in several directions.The flexible modeling framework allows imposing constraints on threshold as well as regression coefficients.In addition, various assumptions about the variance-covariance structure of the errors are supported, by specifying different types of error structures.These include a general correlation, a general covariance, an equicorrelation and an AR(1) error structure.The general error structures can depend on a categorical covariate, while in the equicorrelation and AR(1) structures both numerical and categorical covariates can be employed.Moreover, in addition to the multivariate probit link, we implement a multivariate logit link for the class of multivariate ordinal regression models.This paper is organized as follows: Section 2 provides an overview of the model class and the estimation procedure, including model specification and identifiability issues.Section 3 presents the main functions of the package.A couple of worked examples are given in Section 4. Section 5 concludes.

Model class and estimation
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.
In order to introduce the class of multivariate ordinal regression models considered in this paper, we start with a brief overview on univariate cumulative link models.

Univariate cumulative link models
Cumulative link models are often motivated by the assumption that the observed categories Y i are a categorized version of an underlying latent variable Y i with where β 0 is an intercept term, x i is a p × 1 vector of covariates, β = (β 1 , . . ., β p ) ⊤ is a vector of regression coefficients and ϵ i is a mean zero error term with distribution function F .The link between the observed variable Y i with K categories and the latent variable Y i is given by: where Typical choices for the distribution function F are the normal and the logistic distributions.

Multivariate ordinal regression
Univariate cumulative link models can be extended to a multivariate setting by assuming the existence of several latent variables with a joint error distribution (see, e.g., Varin and Czado 2010;Bhat, Varin, and Ferdous 2010;Kenne Pagui and Canale 2016).Let Y ij denote an ordinal observation and x ij be a p-dimensional vector of covariates for subject i and outcome j, where i = 1, . . ., n and j ∈ 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 sets J and J i , respectively.Following the cumulative link modeling approach, the ordinal response Y ij is assumed to be a coarser version of a latent continuous variable Y ij .The observable categorical outcome Y ij and the unobservable latent variable Y ij are connected by: where r ij is a category out of K j ordered categories and θ j is a vector of suitable threshold parameters for outcome j with the following restriction: 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 Y ij and the vector of covariates x ij : where β j0 is an intercept term, β j = (β j1 , . . ., β jp ) ⊤ is a vector of regression coefficients, both corresponding to outcome j.We further assume the n subjects to be independent.Note that the number of ordered categories K j as well as the threshold parameters θ j and the regression coefficients β j are allowed to vary across outcome dimensions j ∈ J to account for possible heterogeneity across the response variables.
Category-specific regression coefficients.By employing one set of regression coefficients β j for all categories of the jth 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 (McCullagh 1980) and can be relaxed for one or more covariates by allowing the corresponding regression coefficients to be category-specific (see, e.g., Peterson and Harrell 1990).

Link functions.
The dependence among the different responses is accounted for by assuming that, for each subject i, the vector of error terms ϵ i = [ϵ ij ] j∈J i follows a suitable multivariate distribution.We consider two multivariate distributions which correspond to the multivariate probit and logit link functions.For the multivariate probit link, we assume that the errors follow a multivariate normal distribution: ϵ i ∼ N q i (0, Σ i ).A 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 as proposed by O' Brien and Dunson (2004).For a vector z = (z 1 , . . ., z q ) ⊤ , the multivariate logistic distribution function with ν > 2 degrees of freedom, location vector µ and positive-definite dispersion matrix Σ is defined as: where t ν,R is the q-dimensional multivariate t distribution with ν degrees of freedom and correlation matrix R implied by Σ, g ν (x) = t −1 ν (exp(x)/(exp(x) + 1)), with t −1 ν the quantile function of the univariate t distribution with ν degrees of freedom and σ 2 1 , . . ., σ 2 q the diagonal elements of Σ. Hirk, Hornik, and Vana (2019) employed this t copula based multivariate logistic family, while Nooraee, Abegaz, Ormel, Wit, and Van den Heuvel (2016) used a multivariate t distribution with ν = 8 degrees of freedom as an approximation for this multivariate logistic distribution.The employed distribution family differs from the conventional multivariate logistic distributions of Gumbel (1961) or Malik and Abraham (1973) 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.

Identifiability issues
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 Σ i to be a covariance matrix with diagonal elements [σ 2 ij ] j∈J i , only the quantities β j /σ ij and (θ j,r −β j0 )/σ ij are identifiable in the model in Equation 1.Hence, in order to obtain an identifiable model the parameter set is typically constrained in one of the following ways: • Fixing the intercept β j0 (e.g., to zero), using flexible thresholds θ j and fixing σ ij (e.g., to unity) ∀j ∈ J i , ∀i ∈ {1, . . ., n}.
Note that the first two options are the most commonly used in the literature.All of these alternative model parameterizations are supported by the mvord package, allowing the user to choose the most convenient one for each specific application.Table 2 in Section 3.5 gives an overview on the identifiable parameterizations implemented in the package.

Error structures
Different structures on the covariance matrix Σ i can be imposed.

Basic model
The basic multivariate ordinal regression model assumes that the correlation (and possibly variance, depending on the parameterization) parameters in the distribution function of the ϵ i are constant for all subjects i.
Correlation.The dependence between the multiple measurements or outcomes can be captured by different correlation structures.Among them, we concentrate on the following three: • The general correlation structure assumes different correlation parameters between pairs of outcomes COR(ϵ ik , ϵ il ) = ρ kl .This error structure is among the most common in the literature (e.g., Scott and Kanaroglou 2002;Bhat et al. 2010;Kenne Pagui and Canale 2016).
• The equicorrelation structure COR(ϵ ik , ϵ il ) = ρ implies that the correlation between all pairs of outcomes is constant.
• When faced with longitudinal data, especially when moderate to long subject-specific time series are available, an AR(1) autoregressive correlation model of order one can be employed.Given equally spaced time points this AR(1) error structure implies an exponential decay in the correlation with the lag.If k and l are the time points when Y ik and Y il are observed, then COR(ϵ ik , ϵ il ) = ρ |k−l| .
Variance.If a parameterization with identifiable variance is used (see Section 2.3), in the basic model we assume that for each multiple measurement the variance is constant across all subjects (VAR(ϵ ij ) = σ 2 j ).

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.

Correlation.
For each subject i and each pair (k, l) from the set J i , the correlation parameter ρ ikl is assumed to depend on a vector s i of m subject-specific covariates.In this paper we use the hyperbolic tangent transformation to reparameterize the linear term α 0kl + s ⊤ i α kl in terms of a correlation parameter: .
If α kl = 0 for all k, l ∈ J i , this model would correspond to the general correlation structure in the basic model.Moreover, if α 0kl = 0 and α kl = 0 for all k, l ∈ 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: Variance.Similarly, one can model the heterogeneity among the subjects through the variance parameters VAR(ϵ ij ) = σ 2 ij by employing the following linear model on the log-variance: Note that other suitable link functions for the correlation and variance parameterizations could also be applied.The positive-semi-definiteness of the correlation (or covariance) matrix Σ i can be ensured by the use of special algorithms such as the one proposed by Higham (1988).

Composite likelihood 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 (Varin, Reid, and Firth 2011).
For a given parameter vector δ, which contains the threshold parameters, the regression coefficients and the parameters of the error structure, the likelihood is given by: , ) is a Cartesian product, w i are subject-specific non-negative weights (which are set to one in the default case) and f i,q i is the q i -dimensional density of the error terms ϵ 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: Denoting by f i,1 and f i,2 the uni-and bivariate density functions corresponding to the error distribution, the uni-and bivariate probabilities are given by: The maximum pairwise likelihood estimates δpℓ are obtained by direct maximization of the composite likelihood given in Equation 3. 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 Pinheiro and Bates (1996).This parameterization has the advantage that it can be easily applied to correlation matrices.Thirdly, for 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 → ∞.In addition, the maximum pairwise likelihood estimator is asymptotically normal with asymptotic mean δ and a covariance matrix which equals the inverse of the Godambe information matrix: where H(δ) is the Hessian (sensitivity matrix) and V (δ) the variability matrix.The variability matrix V (δ) and the Hessian H(δ) can be estimated as: , and , where pℓ i (δ) is the component of the pairwise log-likelihood corresponding to subject i and pℓ ikl (δ) corresponds to subject i and pair (k, l).
In order to compare different models, the composite likelihood information criterion by Varin and Vidoni (2005) can be used: CLIC where k = 2 corresponds to CLAIC and k = log(n) corresponds to CLBIC).A comprehensive overview and further details on the properties of the maximum composite likelihood estimates are provided in Varin (2008).

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 2.3).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.In this section we assume for sake of notational simplicity that Σ i is a correlation matrix implying that marginally the errors of subject i have variance one and univariate marginal distribution function F 1 for each outcome j.In the more general case with non-constant variances σ 2 ij , F j i,1 should be used instead of F 1 .The marginal cumulative probabilities implied by the model in Equation 1are then given by the following relationship: 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: where f 1 is the density corresponding to F 1 , x ijv is the vth element in x ij and β jv is the vth element in β 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: where all elements of xij are equal to x ij except for the vth element, which is equal to xijv = x ijv + ∆x ijv for the change ∆x ijv in the variable x v .We refer to Greene and Hensher (2010) and Boes and Winkelmann (2006) 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: An increase of one unit in variable v of outcome j (given that all other variables are held constant) changes the probit of the probability that category r ij or lower is observed by the value of the coefficient β jv of this variable.In other words P(Y ij ≤ r ij |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 (McCullagh 1980): For a one unit increase in variable v of outcome j holding all the others constant, we expect a change of size of the coefficient β 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: This means for a one unit increase in variable v of outcome j, holding all the other variables constant, changes the odds ratio by exp(β jv ).In other words, the odds after a one unit change in variable v of outcome j are the odds before the change multiplied by exp(−β jv ): 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 the concept of importance.Reusens and Croux (2017) extend an approach for comparing coefficients of probit and logit models by Hoetker (2007) in order to compare the coefficients across repeated measurements.They analyze the importance ratio where β j,base is the coefficient of a 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 compensation variation: 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.It is to be noted that the importance ratio R jv depends on the scale of the base variable and variable v of outcome j.This implies that the comparison among the measurements j should be done only if the scales of these variables are equal across the multiple measurements.For this purpose, standardization of the covariates for each measurement should be employed.

Implementation
The mvord package contains six data sets and the built-in functions presented in Table 1.
Multivariate ordinal regression models in the R package mvord can be fitted using the main function mvord().Two different data structures can be passed on to the mvord() function through the use of two different multiple measurement objects MMO and MMO2 in the lefthand side of the model formula.MMO uses a long data format, which has the advantage that it allows for varying covariates across multiple measurements.This flexibility requires to specify a subject index as well as a multiple measurement index.In contrast to MMO, the multiple measurement object MMO2 has a simplified data structure but is only applicable in settings where the covariates do not vary between the multiple measurements.In this case, the multiple ordinal observations as well as the covariates are stored in different columns of a 'data.frame'.We refer to this data structure as wide data format.For illustration purposes we use a worked example based on a simulated data set consisting of 100 subjects for which two multiple ordinal responses (Y1 and Y2), two continuous covariates (X1 and X2) and two factor covariates (f1 and f2) are available.The ordinal responses each have three categories labeled with 1, 2 and 3.

Implementation MMO
The fitting function mvord() requires two compulsory input arguments, a formula argument and a data argument:

Data structure
In MMO we use a long format for the input of data, where each row contains a subject index i, a multiple measurement index j, an ordinal observation Y and all the covariates (X1 to Xp).This long format data structure is internally transformed to an n × q matrix of responses which contains NA in the case of missing entries and a list of covariate matrices X j for all j ∈ J.This is performed by the multiple measurement object MMO(Y, i, j) which specifies the column names of the subject index and the multiple measurement index in data.The column containing the ordinal observations can contain integer or character values or inherit from class (ordered) 'factor'.When using the long data structure, this column is basically a concatenated vector of each of the multiple ordinal responses.Internally, this vector is then split according to the measurement index.Then the ordinal variable corresponding to each measurement index is transformed into an ordered 'factor'.For an integer or a character vector the natural ordering is used (ascending, or alphabetical).If for character vectors the alphabetical order does not correspond to the ordering of the categories, the optional argument response.levelsallows to specify the levels for each response 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.If all the multiple measurements use the same number of classes and same labeling of the classes, the column Y can be stored as an ordered 'factor' (as it is often the case in longitudinal studies).
The order of the multiple measurements is needed when specifying constraints on the threshold or regression parameters (see Sections 3.5 and 3.6).This order is based on the type of the multiple measurement index column in data.For 'integer', 'character' or 'factor' the natural ordering is used (ascending, or alphabetical).If a different order of the multiple responses is desired, the multiple measurement index column should be an ordered factor with a corresponding ordering of the levels.

Formula
The multiple measurement object MMO including the ordinal responses Y, the subject index i and the multiple measurement index j is passed on the left-hand side of a 'formula' object.The covariates X1, . . ., Xp are passed on the right-hand side.In order to ensure identifiability intercepts can be included or excluded in the model depending on the chosen model parameterization.

Model without intercept.
If the intercept should be removed, the formula can be specified in the following ways:

. + Xp
Note on the intercept in the formula.We differ in our approach of specifying the model formula from the model formula specification in, e.g., MASS::polr() or ordinal::clm(), in that we allow the user to specify models without an intercept.This option is not supported in the MASS and ordinal packages, where an intercept is always specified in the model 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 2.3.
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 'formula' object and hence set to zero, the rest of the covariates should be specified in such a way that multicollinearity does not arise.This is of primary importance when including categorical covariates, where one category will be taken as baseline by default.

Implementation MMO2
We use the same worked example as above to show the usage of mvord() with the multiple measurement object MMO2.The data set data_mvord_toy has already the required data structure with each response and all the covariates in separate columns.The multiple measurement object MMO2 combines the different response columns on the left-hand side of the formula object: The multiple measurement object MMO2 is only applicable for settings where the covariates do not vary between the multiple measurements.

Data structure
The data structure applied by MMO2 is slightly simplified, where the multiple ordinal observations as well as the covariates are stored as columns in a 'data.frame'.Each subject i corresponds to one row of the data frame, where all outcomes Y i1 , . . ., Y iq (with missing observations set to NA) and all the covariates x i1 , . . ., x ip are stored in different columns.Ideally each outcome column is of type ordered 'factor'.If columns of the responses have types like 'integer', 'character' or 'factor' a warning is displayed and the natural ordering is used (ascending, or alphabetical).

. + Xp
The ordering of the responses is given by the ordering in the left-hand side of the model formula.MMO2 performs like cbind() in fitting multivariate models in, e.g., lm() or glm().

Link functions
The multivariate link functions are specified as objects of class 'mvlink', which 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 are passed on, these will be used in the computation of the standard errors.This design was inspired by the design of the 'family' class in package stats and 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 bivariate normal probabilities which enter the pairwise log-likelihood are computed with package pbivnorm (Genz and Kenkel 2015).The multivariate probit link is the default link function and can be specified by: link = mvprobit() For the multivariate logit link a t copula based multivariate distribution with logistic margins is used (as explained in Section 2.2) and can be specified by: link = mvlogit(df = 8L) The mvlogit() function has an optional integer valued argument 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 ν ≈ 8, the multivariate logistic distribution in Equation 2is well approximated by a multivariate t distribution (O' Brien and Dunson 2004).This is also the value chosen by Nooraee et al. (2016) 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 Fortran code from Alan Genz (Genz and Bretz 2009) 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 df to the user, who could then estimate a separate model for each value in the grid.The best model can be chosen by CLAIC or CLBIC.

Error structures
Different error structures are implemented in mvord and can be specified through the argument error.structure.The error structure objects are of class 'error_struct'.This approach slightly differs from the approach in package nlme, where the error structure is defined by two classes: 'varFunc' for the variance function and 'corStruct' for the correlation structure.We also define the following subclasses for the error structures: 'cor_general' (similar to nlme's 'corSymm'), 'cor_equi' (similar to 'corCompSymm'), 'cor_ar1' (similar to 'corAR1') and 'cov_general' (similar to 'corSymm' with variance function 'varIdent').The different error structures are chosen through the argument error.structure.

Basic model
In the basic model we support three different correlation structures and one covariance structure.
Correlation.For the basic model specification the following correlation structures are implemented in mvord: • cor_general(formula = ~1) -A general error structure, where the correlation matrix of the error terms is unrestricted and constant across all subjects: COR(ϵ ik , ϵ il ) = ρ kl .
Variance.A model with variance parameters VAR(ϵ ij ) = σ 2 j corresponding to each outcome, when the identifiability requirements are fulfilled, can be specified in the following way: • cov_general(formula = ~1) -The estimation of σ 2 j is only implemented in combination with the general correlation structure.

Extending the basic model
The basic model can be extended by allowing covariate dependent error structures.
Correlation.The following correlation structures are implemented in mvord: • cor_general(formula = ~f) -For the heterogeneous general correlation structure, the current implementation only allows the use of one 'factor' variable 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.

Variance.
The following variance structure is implemented in mvord: • cov_general(formula = ~f) -As in the basic model, the estimation of the heterogeneous variance parameters can be performed for the general covariance structure.A subject-specific 'factor' f can be used as a covariate in the log-variance equation.In addition to the correlation matrices, which are estimated for each factor level of f, a vector of dimension q of variance parameters will be estimated for each factor level.

Constraints on 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 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 some of the thresholds have to be fixed for some of the parameterizations presented in Table 2 in order to ensure identifiability of the model.

Threshold constraints across responses
Such constraints can be imposed by a vector of positive integers 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 the worked example, if one wishes to restrict the threshold parameters of the two outcomes Y1 and Y2 to be equal (θ 1 = θ 2 ), this can be specified as: where the first value corresponds to the first response Y1 and the second to the second response Y2.This order of the responses is defined as explained in Sections 3.1 and 3.2

Fixing threshold values
Values for the threshold parameters can be specified by the argument threshold.values.
For this purpose the user can pass a 'list' with q elements, where each element is a 'vector' of length K j − 1 (where K j is the number of ordered categories for ordinal outcome j).A numeric value in this vector fixes the corresponding threshold parameter to a specified value while NA leaves the parameter flexible and indicates it should be estimated.
After specifying the error structure (through the error.structureargument) and the inclusion/exclusion of an intercept in the formula argument, the user can choose among five possible options for fixing the thresholds: • Leaving all thresholds flexible.
• Fixing the first threshold θ j,1 to a constant a j for all j ∈ J.
• Fixing the first and second thresholds θ j,1 = a j , θ j,2 = b j for all outcomes with K j > 2.
• Fixing the first and last thresholds θ j,1 = a j , θ j,K j −1 = b j for all outcomes with K j > 2.
• An extra option is fixing all of the threshold parameters, for all j ∈ J.
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 outcome).Table 2 provides information about the options available for each combination of error structure and intercept, as well as about the default values in case the user does not specify any threshold values.In the presence of binary observations (K j = 2), if a 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, if the following restrictions should apply to the worked example: this can be specified as: threshold.values= list(c(-1, NA), c(-1, NA))

Error structure Intercept
Threshold parameters All flexible One fixed Two fixed Two fixed All fixed Table 2: This table displays different model parameterizations in the presence of ordinal observations (K j > 2 ∀j ∈ J).The row cor includes error structures cor_general, cor_equi and cor_ar1, while row cov includes the error structure cov_general.The minimal restrictions (default) to ensure identifiability are given in green.The default threshold values (in case threshold.values= NULL) are always a j = 0 and b j = 1.

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, 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 the 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 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.

Coefficient constraints across responses
Such constraints can be specified by the argument coef.constraints, which can be either a vector or a matrix of integer values.If vector constraints of the type β k = β l are desired, which should hold for all 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.
Consider the following specification of the latent processes in the worked example: where the regression coefficients for variables X1 and X2 are set to be equal across the two outcomes (β 1 = β 2 ) by: where β jk,r denotes the regression coefficient of covariate k in the linear predictor of the rth cumulative probit or logit for measurement j.By the first restriction, for the first outcome two regression coefficients are employed for covariate X1: β 11,1 for the first linear predictor and β 11,2 for the second linear predictor.Covariate X2 only appears in the model for the second outcome.For each of the two linear predictors a different regression coefficient is estimated: β 22,1 and β 22,2 .
The constraints are set up as a named list where the names correspond to the names of all covariates in the model matrix.To check the name of the covariates in the model matrix one can use the auxiliary function names_constraints() available in mvord (see also next subsection): The number of rows is equal to the total number of linear predictors j (K j − 1) of the ordered responses; in the example above 2 + 2 = 4 rows.The number of columns represents the number of parameters to be estimated for each covariate: coef.constraints = list( X1 = cbind(c(1, 0, 0, 0), c(0, 1, 0, 0), c(0, 0, 1, 1)), X2 = cbind(c(0, 0, 1, 0), c(0, 0, 0, 1)), f2c2 = cbind(rep(1, 4))) For more details we refer the reader to the documentation of the VGAM package.

Interaction terms and categorical covariates
When constraints on the regression coefficients should be specified in models with interaction terms or categorical covariates, the coef.constraintsmatrix 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 names_constraints() before constructing the coef.constraintsand coef.valuesmatrices.The names of each column in the covariate matrix can be obtained by:

R> formula <-MMO2(Y1, Y2) ~1 + X1 : X2 + f1 + f2 * X1 R> names_constraints(formula, data = data_mvord_toy)
[1] "(Intercept)" "f1B" "f1C" "f2c2" [5] "X1" "X1:X2" "f2c2:X1" This should be used when setting up the coefficient constraints.Please note that by default category A for factor f1 and category c1 for factor f2 are taken as baseline categories.This can be changed by using the optional argument contrasts.In models without intercept, the estimated threshold parameters relate to the baseline category and the coefficients of the remaining factor levels can be interpreted as a shift of the thresholds.

weights.name
Weights on each subject i are chosen in a way that they are constant across multiple measurements.Weights should be stored in a column of data.The column name of the weights in data should be passed as a character string to this argument by weights.name= "weights".If weights.name= NULL all weights are set to one by default.Negative weights are not allowed.

offset
If offsets are not specified in the model formula, a list with a vector of offsets for each multiple measurement can be passed.
contrasts can be specified by a named list as in the argument contrasts.arg of the default method of model.matrix().

PL.lag
In longitudinal studies, where q i is possibly large, the pairwise likelihood estimation can be time consuming as it is built from all two-dimensional combinations of j, k ∈ 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 Varin and Czado (2010).Assuming that, for each subject i, we have a time series of consecutive ordinal observations, the ith component of the pairwise likelihood has the following form: The lag can be fixed by a positive integer argument PL.lag and it can only be used along with 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.

Control function mvord.control()
Control arguments can be passed by the argument control and are hidden in the sub-function mvord.control()with the following arguments.

solver
This argument can either be a character string or a function.All general purpose optimizers of the R package optimx (Nash and Varadhan 2011;Nash 2014) can be used for maximization of the composite log-likelihood by passing the name of the solver as a character string to the solver argument.The available solvers in optimx are, at the time of writing, "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "nlm", "nlminb", "spg", "ucminf", "newuoa", "bobyqa", "nmkb", "hjkb", "Rcgmin" and "Rvmmin".The default in mvord is solver "newuoa".The "BFGS" solver performs well in terms of computation time, but it suffers from convergence problems, especially for the mvlogit() link.
Alternatively, the user has the possibility of applying other solvers by using a wrapper function with arguments starting.valuesand objFun of the following form: solver = function(starting.values,objFun) { optRes <-solver.function(...) list(optpar = optRes$optpar, objvalue = optRes$objvalue, convcode = optRes$convcode, message = optRes$message) } solver.function()should return a list with the following three elements: optpar, objvalue and convcode.The element optpar should be a vector of length equal to the number of parameters to optimize containing the estimated parameters, while the element objvalue should contain the value of the objective function after the optimization procedure.The convergence status of the optimization procedure should be returned in element convcode with 0 indicating successful convergence.Moreover, an optional solver message can be returned in element message.

solver.optimx.control
A list of control arguments that are to be passed to the function optimx().For further details see Nash and Varadhan (2011).

se
If se = TRUE standard errors are computed analytically using the Godambe information matrix (see Section 2.5).

start.values
A list of starting values for threshold as well as regression coefficients can be passed by the argument start.values.This list contains a list (with a vector of starting values for each dimension) theta of all flexible threshold parameters and a list beta of all flexible regression parameters.

Output and methods for class 'mvord'
The function mvord() returns an object of class 'mvord', which is a list containing the following components: • beta: A named 'matrix' of regression coefficients.
• theta: A named 'list' of threshold parameters.
• sebeta: A named 'matrix' of standard errors of the regression coefficients.
• setheta: A named 'list' of standard errors of the threshold parameters.
• seerror.struct:A 'vector' of standard errors for the parameters of the error structure.
• rho: A 'list' of objects that are used in mvord().
Several methods are implemented for the class 'mvord'.These methods include a summary() and a print() function to display the estimation results, a coef() function to extract the regression coefficients, a thresholds() function to extract the threshold coefficients and a function error_structure() to extract the estimated parameters of the correlation/covariance structure of the errors.The pairwise log-likelihood can be extracted by the function logLik(), function vcov() extracts the variance-covariance matrix of the parameters and nobs() provides the number of subjects.Other standard methods such as terms() and model.matrix()are also available.Functions AIC() and BIC() can be used to extract the composite likelihood information criteria CLAIC and CLBIC.

Examples
In credit risk applications, multiple credit ratings from different credit rating agencies are available for a panel of firms.Such a data set has been analyzed in Hirk et al. (2019), where a multivariate model of corporate credit ratings has been proposed.Unfortunately, this original data set is not freely re-distributable.Therefore, we resort to the simulation of illustrative data sets by taking into consideration key features of the original data.
We simulate relevant covariates corresponding to firm-level and market financial ratios in the original data set.The following covariates are chosen in line with literature on determinants of credit ratings (e.g., Campbell, Hilscher, and Szilagyi 2008;Puccia, Collett, Kernan, Palmer, Mettrick, and Deslondes 2013): LR (liquidity ratio relating the current assets to current liabilities), LEV (leverage ratio relating debt to earnings before interest and taxes), PR (profitability ratio of retained earnings to assets), RSIZE (logarithm of the relative size of the company in the market), BETA (a measure of systematic risk).We fit a distribution to each covariate using the function fitdistr() of the MASS package.The best fitting distribution among all available distributions in fitdistr() has been chosen by AIC.
We generate two data sets for illustration purposes.The first data set data_cr consists of multiple ordinal rating observations at the same point in time for a collection of 690 firms.We generate ratings from four rating sources rater1, rater2, rater3 and rater4.Raters rater1 and rater2 assign ordinal ratings on a five-point scale (from best to worst A, B, C, D and E), rater3 on a six-point scale (from best to worst, F, G, H, I, J and K) and rater4 distinguishes between investment and speculative grade firms (from best to worst, L and M).The panel of ratings in the original data set is 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 set.For rater1 we remove 5%, for rater2 30%, and for rater3 50% of the observations.This data set has a wide data format.
The second data set data_cr_panel contains ordinal rating observations assigned by one rater to a panel of 1415 firms over a period of eight years on a yearly basis.In addition to the covariates described above, a business sector variable (BSEC) with eight levels is included for each firm.This data set has a long format, with 11320 firm-year observations.

Example 1: A simple model of firm ratings assigned by multiple raters
The first example presents a multivariate ordinal regression model with probit link and a general correlation error structure cor_general(~1).The simulated data set contains the ratings assigned by raters rater1, rater2, rater3 and rater4 and the five covariates LR, LEV, PR, RSIZE and BETA for a cross-section of 690 firms.A value of NA indicates a missing observation in the corresponding outcome variable.R> data("data_cr", package = "mvord") R> head(data_cr, n = 3)  RSIZE BETA 1 -6.365053 0.8358773 2 -7.839813 0.4895358 3 -7.9766500.8022900 R> str (data_cr, vec.len = 2.9) 'data.frame':690 obs. of 10 variables: $ rater1 : Ord.factor w/ 5 levels "A"<"B"<"C"<"D"<..: 2 3 3 2 5 4 3 ... $ rater2 : Ord.factor w/ 5 levels "A"<"B"<"C"<"D"<..: 2 4 4 2 5 NA 3 ... $ rater3 : Ord.factor w/ 6 levels "F"<"G"<"H"<"I"<..: 3 NA NA NA 6 NA 2 ... $ rater4 : Ord.factor w/ 2 levels "L"<"M": 1 2 2 1 2 2 2 ... $ firm_id: int 1 2 3 4 5 6 7  The threshold parameters are labeled with the name of the corresponding outcome and the two adjacent categories which are separated by a vertical bar |.For each covariate the estimated coefficients are labeled with the covariate name and a number.This number is from the sequence along the number of columns in the list element of constraints() which corresponds to the covariate.Note that if no constraints are set on the regression coefficients, this number of the coefficient corresponds to the outcome dimension.If constraints are set on the parameter space, we refer the reader to Section 4.2.The last part of the summary contains the estimated error structure parameters.For error structures cor_general and cov_general the correlations (and variances) are displayed.The coefficients corresponding to the error structure are displayed for cor_ar1 and cor_equi.Correlations and Fisher-z scores for each subject are obtained by function error_structure().

Example 2: A more elaborate model of ratings by multiple raters
In the second example, we extend the setting of Example 1 by imposing constraints on the regression as well as on the threshold parameters and changing the link function to the multivariate logit link.We include the following features in the model: • We assume that rater1 and rater2 use the same rating methodology.This means that they use the same rating classes with the same labeling and the same thresholds on the latent scale.Hence, we set the following constraints on the threshold parameters: threshold.constraints= c(1, 1, 2, 3) • We assume that some covariates are equal for some of the raters.We assume that the coefficients of LR and PR are equal for all four raters, that the coefficients of RSIZE are equal for the raters rater1, rater2 and rater3 and the coefficients of BETA are the same for the raters rater1 and rater2.The coefficients of LEV are assumed to vary for all four raters.These restrictions are imposed by: The estimation can now be performed by the function mvord(): R> res_cor_logit <-mvord(formula = MMO2(rater1, rater2, rater3, rater4) + 0 + LR + LEV + PR + RSIZE + BETA, data = data_cr, link = mvlogit(), + coef.constraints= cbind(LR = c(1, 1, 1, 1), LEV = c(1,2,3,4), The results are displayed by the function summary():

R> summary(res_cor_logit, call = FALSE)
Formula: MMO2(rater1, rater2, rater3, rater4 If constraints on the threshold or regression coefficients are imposed, duplicated estimates are not displayed.If thresholds are set equal for two outcome dimensions only the thresholds for the former dimension are shown.In the example above only the thresholds for rater1 are displayed.For each covariate the estimated coefficients are labeled with the covariate name and a number.This number is from a sequence along the number of columns in the list element of the corresponding covariate in constraints() (see Section 3.6).The auxiliary function constraints() can be used to extract the constraints on the coefficients.The column names of the constraint matrices for each outcome correspond to the coefficient names displayed in the summary.For each covariate the coefficients to be estimated are numbered consecutively.In the above example this means that for covariates LR and PR only one covariate is estimated, a coefficient for each outcome is estimated for LEV, while for covariate RSIZE two and for covariate BETA three coefficients are estimated.For example, the coefficient BETA 1 is used by rater1 and rater2, the coefficient BETA 2 is used by rater3 while BETA 3 is the coefficient for rater4.The constraints for covariate BETA can be extracted by:

Comparing the model fits of examples one and two
Note that the composite likelihood information criteria can be used for model comparison.For objects of class 'mvord' CLAIC and CLBIC are computed by AIC() and BIC(), respectively.The value of the pairwise log-likelihood of the two models can be extracted by logLik().
We include five financial ratios as covariates in the model with an intercept by a formula with multiple measurement object MMO: formula = MMO(rating, firm_id, year) ~LR + LEV + PR + RSIZE + BETA Additionally, the model has the following features: • The threshold parameters are constant over the years.This can be specified through the argument threshold.constraints: threshold.constraints= rep(1, nlevels(data_cr_panel$year)) • In order to ensure identifiability in a model with intercept, some thresholds need to be fixed.We fix the first thresholds for all outcome dimensions to zero by the argument threshold.values: threshold.values= rep(list(c(0, NA, NA, NA)), 8) • We assume that there is a break-point in the regression coefficients after year4 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 year1, year2, year3 and year4 and a different set for year5, year6, year7 and year8.This can be specified through the argument coef.constraints: coef.constraints = c(rep(1, 4), rep(2, 4)) • Given the longitudinal aspect of the data, an AR(1) correlation structure is an appropriate choice.Moreover, we use the business sector as a covariate in the 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.

Conclusion
The present paper is meant to provide a general overview on the R package mvord, which implements the estimation of multivariate ordinal probit and logit regression models using the pairwise likelihood approach.We offer the following features which (to the best of our knowledge) enhance the currently available software for multivariate ordinal regression models in R: • Different error structures like a general correlation and a covariance structure, an equicorrelation structure and an AR(1) structure are available.
• We 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.
• We allow for outcome-specific threshold parameters.
• We allow for outcome-specific regression parameters.
• 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 outcomes).
• We offer the possibility to choose different parameterizations, which are needed in ordinal models to ensure identifiability.
Additional flexibility is achieved by allowing the user to implement alternative multivariate link functions or error structures (e.g., alternative transformations for the variance or correlation parameters can be implemented).Furthermore, the long as well as the wide data format are supported by either applying MMO or MMO2 as a multiple measurement object to estimate the model parameters.The functionality of the package is illustrated by a credit risk application.Further examples from different areas of application are presented in the package vignette.
Further research and possible extensions of mvord could consist of the implementation of variable selection procedures in multivariate ordinal regression models and the inclusion of multivariate semi-or non-parametric ordinal models.

Figure 1 :
Figure 1: This figure displays the rating distribution of all the raters.

Figure 2 :
Figure2: This figure displays agreement plots of the predicted categories of the model res_cor_logit against the observed rating categories for all raters.For each observed rating class the distribution of the predicted ratings is displayed.

For
the fixed threshold coefficient year1 A|B, the z values and the corresponding p values are set to NA.The default error_structure() method for a 'cor_ar1' gives: -0.55231725 -0.06202790 -0.08559155 -0.06659849 V21 V22 V23 -0.68342920 -0.86391149 -0.75799729In addition, the correlation parameters ρ i for each firm are obtained by choosing type = "corr" in error_structure():R> head(error_structure(res_AR1_probit, type = "corr"), n = 3) correlation matrices for each specific firm are obtained by choosing type = "sigmas" in error_structure(): (McCullagh 1980tz 2012rs on the latent scale (see, e.g.,Agresti 2010;Tutz 2012).In such a setting the ordinal response variable Y i follows a multinomial distribution with parameter π i .Let denote by π ir i the probability that observation i falls in category r i .Then the cumulative link model(McCullagh 1980) is specified by: The threshold coefficients can be extracted by the function thresholds():

Table 3 :
The model fit of examples one and two are compared by means of BIC and AIC.From Table 3 we observe that the model of Example 2 has a lower BIC and AIC, which indicates a better model fit.This table displays measures of fit for the multivariate probit model in Example 1 (presented in Section 4.1) and the multivariate logit model in Example 2 (presented in Section 4.2).