Regularized Ordinal Regression and the ordinalNet R Package

Regularization techniques such as the lasso (Tibshirani 1996) and elastic net (Zou and Hastie 2005) can be used to improve regression model coefficient estimation and prediction accuracy, as well as to perform variable selection. Ordinal regression models are widely used in applications where the use of regularization could be beneficial; however, these models are not included in many popular software packages for regularized regression. We propose a coordinate descent algorithm to fit a broad class of ordinal regression models with an elastic net penalty. Furthermore, we demonstrate that each model in this class generalizes to a more flexible form, for instance to accommodate unordered categorical data. We introduce an elastic net penalty class that applies to both model forms. Additionally, this penalty can be used to shrink a non-ordinal model toward its ordinal counterpart. Finally, we introduce the R package ordinalNet, which implements the algorithm for this model class.


Introduction
Ordinal regression models arise in contexts where the response variable belongs to one of several ordered categories (such as 1="poor", 2="fair", 3="good", 4="excellent").One of the most common regression models for this type of data is the cumulative logit model (McCullagh, 1980), which is also known as the proportional odds model or the ordinal logistic regression model.Other ordinal regression models include the stopping ratio model, the continuation ratio model, and the adjacent category model.The VGAM R package (Yee and Wild, 1996;Yee, 2010Yee, , 2015) ) fits all of the aforementioned models, but without regularization or variable selection.The SAS CATMOD procedure also fits some of these models (SAS Institute Inc, 2017).Popular CRAN packages for penalized regression, such as penalized (Goeman et al., 2014) and glmnet (Friedman et al., 2010), do not currently fit ordinal models.
Some algorithms and software already exist for penalized ordinal regression models.The R package lrm (Harrell, 2015) fits the cumulative logit model with quadratic (ridge regression) penalty.The R packages glmnetcr (Archer, 2014a) and glmpathcr (Archer, 2014b) fit stopping ratio models with the elastic net penalty.Archer (2014a,b) refers to these as continuation ratio models, but we define stopping ratio and continuation ratio models in the same way as Yee (2010).Archer et al. (2014) also implemented the generalized monotone incremental forward stagewise (GMIFS) algorithm for regularized ordinal regression models in the R package ordinalgmifs.This procedure finds a solution path similar to the L 1 norm (lasso) penalty.In fact it is the same solution path if the lasso path is monotone for each coefficient, but in other cases the GMIFS and lasso solution paths differ (Hastie et al., 2009).Some drawbacks of this algorithm are that it fits a single solution path and does not have the flexibility of the elastic net mixing parameter (usually denoted by α).It can also be computationally expensive because the entire solution path must be fit in small increments, whereas the lasso and elastic net solution path can be obtained only at specified values of the regularization tuning parameter (usually denoted by λ).A sequence of, say, twenty values may be enough to tune a model by cross validation and will usually be faster than fitting a longer sequence.
To summarize, algorithms for ordinal regression either do not allow regularization, or they apply to specific models.Hence, options are limited for ordinal regression with a large number of predictors.In that context, our contribution to this growing body of software and literature is threefold.First, we propose a general coordinate descent algorithm to fit a rich class of multinomial regression models, both ordinal and non-ordinal, with elastic net penalty.
Second, we define a class of models that (a) can be fit with the elastic net penalty by the aforementioned algorithm, (b) contains some of the most common ordinal regression models, (c) is convenient for modularizing the fitting algorithm, and (d) has both a parallel (ordinal) and a nonparallel form for each model (discussed in the next paragraph).We call this the elementwise link multinomial-ordinal (ELMO) class of models.This class is a subset of vector generalized linear models (Yee, 2015).Each model in this class uses a multivariate link function to link multinomial probabilities to a set of linear predictors.The link function can be conveniently written as a composite of two functions.The first determines the model family (e.g.cumulative probability, stopping ratio, continuation ratio, or adjacent category).The second is a standard link function (e.g.logit, probit, or complementary log-log), which is applied elementwise to the result of the first function.
Another feature of the ELMO class is that each model has a form that is appropriate for ordinal response data, as well as a more flexible form that can be applied to either ordinal or unordered categorical responses.We will refer to the first as the parallel form and the second as the nonparallel form.For the parallel form, the linear predictors of a given observation only differ by their intercept values-the other coefficients are the same.This restriction is what Yee (2010) refers to as the parallelism assumption.The nonparallel form allows all of the coefficients to vary.An example from the ELMO class is the proportional odds model, which is a parallel model that has a nonparallel counterpart, the partial proportional odds model (Peterson et al., 1990).For more details, see Section 2.5.
Finally, we propose an elastic net penalty class that applies to both the parallel and nonparallel forms.It can also be used to shrink the nonparallel model toward its parallel counterpart.This can be useful in a situation where one would like to fit an ordinal model but relax the parallelism assumption.This can be achieved by over-parameterizing the nonparallel model to include both the nonparallel and parallel coefficients.We call this alternate parametrization the semi-parallel model.Although the regression model itself is not identifiable under this parametrization, the penalized likelihood has a unique optimum (or almost unique, as discussed in Appendix A).
We provide an outline for the remainder of the work.Section 2 defines the ELMO class with specific examples.We also define the parallel, nonparallel, and semi-parallel parametrizations with the elastic net penalty.Section 3 provides the proposed algorithm for fitting multinomial regression models with the elastic net penalty.Section 4 presents a simulation study to compare prediction accuracy of the penalized parallel, nonparallel, and semi-parallel models.Section 5 demonstrates the use of penalized ELMO class models for out-of-sample prediction and variable selection alongside other methods.Section 6 provides details about the ordinalNet R package, which is available on the Comprehensive R Archive Network.Section 7 provides a demonstration of the ordinalNet package.Section 8 contains a summary of the findings and contribution.

Elementwise link multinomial-ordinal (ELMO) class
This section is organized as follows.Section 2.1 introduces commonly used notation.Section 2.2 is the heart of Section 2, defining the ELMO model class.The remaining subsections then discuss particular elements of the ELMO class and issues related to elastic net penalization of this class.Sections 2.3 and 2.4 provide details for the family function and elementwise link function, respectively.The parallel, nonparallel, and semi-parallel forms are discussed in Sections 2.5 and 2.6.Section 2.7 discusses the elastic net penalty and formulates the objective function under the three model forms.Finally, Section 2.8 describes the Jacobian of the inverse link function for the ELMO class.

Notation
We introduce commonly used notation.Paper-specific notation is developed throughout the work.For a vector x, we use x T to denote its transpose.For a matrix A we use bracket notation to indicate elements, so that [A] ij , indicates the (i, j) th element of A. ½ K = ½ denotes the length-K column vector of ones and I K×K = I denotes the K ×K identity matrix.In both cases, the dependence on K will be suppressed when it is clear from the context.We use ∇ for the gradient operator and D for the Jacobian operator.Consider a vector-valued function f with vector argument x.As is standard, the Jacobian of f is defined as

An Introduction to the ELMO model class
We now define the ELMO model class.Models within this class are completely specified by their multivariate link function, which is a composite of two functions.The first function determines the model family (e.g.cumulative probability, stopping ratio, continuation ratio, or adjacent category).We refer to these as multinomial-ordinal (MO) families because each has a parallel form specifically for ordinal data, as well as a nonparallel form for any multinomial data, ordinal or unordered.The second function is an elementwise link (EL) function, which applies a standard link function on (0, 1) → R (e.g.logit, probit, or complementary log-log) elementwise to the result of the first function.McCullagh (1980), Wilhelm et al. (1998), and Yee (2010) provide more background on categorical regression models. Let ∼ Multinomial(n i , p i1 , p i2 , . . ., p i(K+1) ) for i = 1, . . ., N .Observations (e.g.subjects or patients) are indexed by i, and N is the total number of observations.Here, x i is an observed vector of covariates (without an intercept), and y i = (y i1 , . . ., y i(K+1) ) T is a random vector of counts summing to n i .The conditional distribution represents n i independent trials which fall into K + 1 classes with probabilities (p i1 , p i2 , . . ., p i(K+1) ) that are a function of x i .The K + 1 probabilities sum to one, so they can be parametrized by the vector p i = (p i1 , p i2 , . . ., p iK ) T .
Let P be the length of x i and let B be a P × K matrix of regression coefficients.Let b 0 be a vector of K intercept values.The covariates are mapped to a vector of K linear predictors, η i , by the relationship η i = b 0 +B T x i .Class probabilities are linked to the linear predictors by η i = g(p i ), where g : S K → R K is a multivariate invertible link function and S K = {p : p ∈ (0, 1) K , p 1 < 1}.Furthermore, g is a composite of two functions, g EL : (0, 1) K → R K and g MO : S K → (0, 1) K .More specifically, ELMO class models have a link function of the form

Family function
The function g MO determines the family of multinomial-ordinal models, such as cumulative probability, stopping ratio, continuation ratio, or adjacent category.
In order to belong to a multinomial-ordinal family, the function g MO must be invertible and have the following Monotonicity Property.This Monotonicity Property ensures that all parallel models in the ELMO class are in fact ordinal models (discussed in Section 2.5).Examples of MO families are given in Table 1.
Definition (Monotonicity Property) For any p ∈ S K , define γ j (p) for j = 1, . . ., K as the sum of the first j elements (i.e.cumulative probabilities). Define Recall that p i = (p i1 , p i2 , . . ., p iK ) T , and let r(p i ) = (p i(K+1) , p iK , . . ., p i2 ) denote the class probabilities in reverse order, leaving out class 1 instead of class K + 1.If g MO is a MO function with Property 1, then the (g MO • r) is a MO function with Property 2 and vice versa.We can refer to one as the "forward" family and the other as the "backward" family.Although the terms "forward" and "backward" are commonly used in the literature, they do not have a consistent interpretation.We follow the naming conventions used in Yee (2010).By these definitions, the forward cumulative probability and stopping ratio families have Property 1, and the forward continuation ratio and adjacent category families have Property 2.
In addition, if g MO is an MO function, then g * MO (p) = 1 − g MO (p) is also an MO function.For the cumulative probability and adjacent category families, this is simply a transformation between the forward and backward families.On the other hand, applying this transformation to the forward (backward) stopping ratio family yields the forward (backward) continuation ratio family.

Elementwise link function
The elementwise link function g EL * : (0, 1) → R must be a monotone, invertible function.It can be almost any link function used for binary data regression, such as logit, probit, or complementary log-log.An important property that some elmentwise link functions satisfy is symmetry, that is For example, logit and probit are symmetric, but complementary log-log is not.Under symmetry, the following model pairs are equivalent, with reversed signs on the coefficients: 1) cumulative probability forward and backward models; 2) adjacent category forward and backward models; 3) forward stopping ratio and forward continuation ratio; and 4) backward stopping ratio and backward continuation ratio.However, if g EL * is not symmetric, such as the complementary log-log, then none of these equivalences hold.

Parallel and nonparallel forms
Each model in the ELMO class has a parallel form and a nonparallel form.The difference between them is that the parallel form restricts the columns of B to be identical.This restriction is referred to as the parallelism assumption by Yee (2010).Let b denote the common column vector and consider the distribution of y i |x i = x i .The Monotonicity Property of g MO , along with the monotonicity requirement on g EL * , ensures that a change in b T x i will shift all cumulative class probabilities in the same direction.This is the defining characteristic of an ordinal regression model.
In contrast, the nonparallel form places no restriction on B, and it does not force the cumulative class probabilities to "shift together" in any way.The nonparallel form is appropriate for unordered multinomial data, although it can also be used as a more flexible model for ordinal data.
A word of caution: the nonparallel cumulative probability model must a have linear predictor vector, B T x, that is monotone to ensure that the cumulative probabilities are non-decreasing.For example, B T x must be monotone increasing for the forward model and monotone decreasing for the backward model.Thus, B should be constrained such that B T x is monotone for any feasible x in the population of interest.This constraint is difficult to implement in practice, especially because the range of feasible x values may not be known.It is more practical to constrain B T x to be monotone for all x in the training sample.However, this may lead to non-monotone probabilities for out-of-sample x, so it is important to be mindful of this.This is not a concern for the parallel cumulative probability model because the MLE (or penalized MLE) will always have monotone intercepts, and hence monotone linear predictors for all x.The other families in Table 1 do not have any restriction on the parameter space.

Semi-parallel form
In most applications with ordinal response data, domain knowledge does not make it clear whether to use the parallel or nonparallel form.With enough observations, one could simply estimate the nonparallel model by maximum likelihood and obtain a good fit.After all, the nonparallel model includes the parallel model as a special case, and the parallel model will give inconsistent coefficient estimates if it is incorrect.
When the number of predictors is large relative to the number of observations, a regularization method such as lasso or elastic net is required.In this case, it is not possible to estimate each coefficient with a high degree of accuracy-a more realistic modeling goal is to build a model for out-of-sample prediction and determine the most important predictors.In this context, one might forgive some incorrectness of the parallel model if it is accurate enough to accomplish the modeling goals.Even if a nonparallel model were the true data generating mechanism, the regularized parallel model could still outperform the regularized nonparallel model for prediction.
The question becomes: how "parallel" does the data need to be to make the parallel model a better choice?If the response categories have a natural ordering, then it seems prudent to leverage this fact by using an ordinal regression model.However, the fact alone that the response is ordinal does not mean that a parallel regression model will be a good fit.Therefore, it also seems prudent to use a model that is sufficiently flexible to allow deviation from the strict parallelism assumption.
With this motivation, we propose a model that ( 1) is ordinal in nature and (2) allows deviation from the parallelism assumption.We call this the semiparallel model.Recall that the nonparallel model specifies η i = b 0 + B T x i , where b 0 is a vector of K intercepts and B is an unrestricted P × K matrix of coefficients.The parallel model restricts the columns of B to be identical and can be parametrized as It is the nonparallel model but overparametrized to include both the parallel and nonparallel coefficients.With the elastic net penalty, the penalized likelihood has a unique solution in most cases (see Appendix A for details).We use the term semi-parallel because for some covariates, the penalized semi-parallel model fit might only contain the parallel coefficient, with the nonparallel coefficients all set to zero.For other covariates, the fit might contain both parallel and nonparallel coefficients.

Elastic net penalty
This section discusses the elastic net penalty for ELMO models.There are many useful resources on regularization, penalized regression, and variable selection, which provide further details in various settings (Hastie et al., 2009;Bickel and Li, 2006;Hesterberg et al., 2008;Friedman et al., 2010;Schifano et al., 2010;Vidaurre et al., 2013).
If the sample size is large enough, it may be possible to accurately estimate a regression model by maximum likelihood.However, in many applications the sample size is not large enough to obtain reliable or even unique estimates.In situations like this, it may be advantageous to optimize a penalized version of the log likelihood function.One such penalty is the elastic net (Zou and Hastie, 2005), which is a generalization of both the lasso (Tibshirani, 1996) and ridge regression penalties.
Lasso and ridge regression are techniques that minimize a penalized likelihood objective function, defined as the negative log-likelihood plus a penalty term that is a function of the coefficient vector.For lasso, the penalty is proportional to the L 1 norm of the coefficient vector, and for ridge regression it is proportional to the squared L 2 norm.Both of these penalties result in a coefficient estimate that is closer to zero than the maximum likelihood estimator, i.e. the estimate is "shrunk" toward zero.This biases the estimates toward zero, but the trade-off is a reduction in variance which often reduces the overall mean squared error.The lasso also has the property that some of the coefficient estimates are shrunk to zero exactly.This provides a natural way to perform variable selection because only the the predictors most associated with the response variable will have nonzero coefficients.
The elastic net penalty is a weighted average between the lasso and ridge regression penalties, and it shares the lasso property of shrinking some coefficients to zero exactly.The weighting parameter, typically denoted by α, must be selected or tuned on the data set.The degree of penalization is controlled by another tuning parameter, typically denoted by λ.Typical practice is to fit the penalized model for a sequence of λ values and use a tuning procedure to select the best value (Hesterberg et al., 2008;Hastie et al., 2009;Arlot et al., 2010;Sun et al., 2013).One tuning procedure is to select the model with best fit as determined by C p , AIC, BIC, or another fitness measure.An alternative approach is using cross-validation to select the value that gives the best out-of-sample prediction.We use cross-validation.
Let b j be the j th element of b and B jk be the element in the j th row and k th column of B. Let N * = N i=1 n i be the total number of multinomial trials in the data set.Let ℓ(•) denote the log-likelihood function for each ELMO model form.Below, we write the elastic net objective function for each model form.

Parallel model
The objective function is

Nonparallel model
The objective function is

Semi-parallel model
The objective function is Here, λ ≥ 0 and α ∈ [0, 1] are the tuning parameters previously described.Also, ρ ≥ 0 is a tuning parameter that determines the degree to which the parallel terms are penalized.Fixing ρ at a very large value will set all parallel coefficients to zero, which is equivalent to the nonparallel model.Fixing λ at a very large value and choosing ρ such that λρ = λ * is equivalent to the parallel model with regularization parameter λ * .Hence, the semi-parallel model includes both the parallel and nonparallel models as special cases.Fixing ρ = 0 will leave the parallel coefficients unpenalized, so the fit will shrink from the maximum likelihood nonparallel model fit toward the maximum likelihood parallel model fit as λ increases from zero.
We follow the common convention of scaling the negative log-likelihood by the number of observations (Hesterberg et al., 2008).This way, when fitting a model to a sample from a given population, a given λ value will have roughly the same degree of penalization regardless of the sample size.This is convenient when tuning λ by cross validation because the tuning data set may have a different sample size than the training data set.We define the sample size as N * rather than N so the model fit is invariant to whether multinomial trials are grouped into a single observation or split into multiple observations with n i = 1 and identical x.

Inverse link function Jacobian
The Jacobian of the inverse link function is required for the coordinate descent algorithm.This computation can be compartmentalized for link functions in the ELMO class because of their composite form.Define h, h EL , h EL * , and h MO to be the inverses of g, g EL , g EL * , and g MO , respectively.The inverse link function can be written as The Jacobian of the inverse link can be written as where The inverse and its derivative are well-known for common elementwise link functions, so we do not discuss these any further (see, e.g., the make.linkfunction in the R package stats).To calculate h(η), it only remains to calculate h MO (δ) and Dh MO (δ).These calculations are presented for specific MO families in Appendix B.

Coordinate descent optimization algorithm
We propose optimizing ELMO class models with the elastic net penalty using a coordinate descent algorithm.Our algorithm mirrors that of Friedman et al. (2007Friedman et al. ( , 2010) ) for generalized linear models.The algorithm is iterative and has an outer and inner loop.The outer loop constructs a quadratic approximation to the log-likelihood-the same quadratic approximation used for iteratively reweighted least squares (IRLS).This approximation is the second order Taylor expansion at the current coefficient estimates.In the spirit of Fisher scoring, the Hessian matrix is replaced by its expectation, the negative Fisher information matrix.This approximation is used as a replacement for the true likelihood in the elastic net objective function, resulting in an expression that can be optimized by coordinate descent.The inner loop cycles through the coefficient estimates, updating each one with the value that marginally optimizes the approximate objective function.Wilhelm et al. (1998) demonstrated the use of IRLS to obtain the maximum likelihood estimates for a broad class of multinomial regression models.This class includes ELMO models but is even more general.This algorithm can easily be applied to any multinomial regression model that links a vector of K probabilities to a vector of K linear combinations of covariates.The link function g does not need to be a composite function or have the Monotonicity Property of the ELMO class.It simply needs to have an inverse h.To apply the coordinate descent algorithm to another model, all that is required is to derive the Jacobian of h.Although ELMO is a rich class of models, there are multinomial regression models outside this class (e.g.multinomial logistic regression).The coordinate descent algorithm is very general.One could use the basic ideas for fitting an elastic net penalized multinomial regression model that is outside the ELMO class.
The rest of this section is organized as follows.In order to formulate the IRLS quadratic approximation, it is more convenient to parametrize ELMO models with a single coefficient vector instead of b 0 , b, and B. Section 3.1 discusses this alternative parameterization.Section 3.2 discusses the elastic net penalty under the alternative parametrization.In Section 3.3 we derive the general form of the score and information matrix for multinomial regression models.Section 3.4 discusses the outer loop of the optimization procedure, which updates the quadratic approximation to the log-likelihood.Section 3.5 discusses the coordinate descent inner loop.Section 3.6 discusses computational efficiency and numerical stability for the coordinate descent updates.
Sections 3.7, 3.8, and 3.9 discuss different aspects of the algorithm specifications.Specifically, Section 3.7 presents a method for choosing a sequence of regularization parameter values for the solution path.Section 3.8 presents a method for choosing starting coefficient values for the optimization algorithm.And Section 3.9 suggests a stopping rule for terminating the algorithm.Section 3.10 summarizes the algorithm in outline form.Finally, Section 3.11 discusses specific optimization issues that can arise with the cumulative probability model family.

ELMO parametrization with a single coefficient vector
Until now, we have written ELMO coefficients in a compact form using an intercept vector b 0 , a coefficient vector b, and a coefficient matrix B. For coordinate descent, it is more convenient to write the model with a single coefficient vector β.To do this, we need to introduce a covariate matrix X i , which is a function of x i .The vector of K linear combinations, η i , can then be written as η i = X T i β for any of the three ELMO model forms.
Let B j• denote the transpose of the j th row of B. The form of X i and β is given below for the parallel, nonparallel, and semi-parallel models.

Parallel model
Nonparallel model .

Elastic net penalty
Suppose β is length Q and let β j denote the j th element.We write the elastic net objective function as where λ > 0 and 0 ≤ α ≤ 1.For ELMO models, c j = 0 for all intercept terms.For the parallel and nonparallel models, c j = 1 for all non-intercept coefficients.
For the semi-parallel model, c j = ρ for the parallel non-intercepts (b), and c j = 1 for the nonparallel non-intercepts (B).These are not firm rules regarding c j , as there may be situations where one wishes to modify the c j to accommodate more elaborate penalization schemes.For example, one might wish leave to some covariates unpenalized or to penalize them with varying degrees.The only requirement is that the c j be nonnegative.Typically, each covariate is standardized to have its sample standard deviation equal to one so that the scale of a covariate does not affect the degree to which its coefficient is penalized (Hesterberg et al., 2008).However, this is not a requirement.

Score and information matrix
In this section, we derive the general form of the score and information matrix for multinomial regression models.The log-likelihood of an observation with probability vector p i can be written as Note that we drop the multinomial term log ni yi1,yi2,...,y i(K+1) because it does not depend on the unknown coefficients, and hence does not affect the model fit.The log-likelihood as a function of β is The score function can be obtained by a chain rule decomposition: where • ½ , and The Fisher information matrix is given by where Because the y i are independent, the full data log-likelihood, score, and information are defined as

Optimization outer loop (quadratic approximation)
The optimization outer loop updates the quadratic approximation to the loglikelihood using a Taylor expansion around the current coefficient estimates.Let β(r) denote the coefficient estimates after the r th outer loop iteration, and let the (r) superscript denote terms that depend on β(r) .Of course, starting values are required for the first iteration, and this topic will be discussed later.
We define the quadratic approximation of ℓ(β) as where z r) .This is the second order Taylor series expansion of ℓ at β(r) , up to an additive constant that does not depend on β.Also, the Hessian is replaced by its expectation, −I( β(r) ).The derivation is provided in Appendix C. The inner loop computes β(r+1) by optimizing the penalized quadratic approximation.

Optimization inner loop (coordinate descent)
For unpenalized maximum likelihood estimation, the quadratic approximation is maximized by the usual weighted least squares solution β r) .With the elastic net penalty, we can still follow the IRLS procedure, but the optimization step no longer has a closed form solution.This is because partial derivatives of the elastic net penalty do not exist at zero for any of the penalized coefficients.The optimization step can instead be done with a coordinate descent procedure.This involves cycling through the coefficient estimates, updating each one with the value that marginally optimizes the approximate objective function.The cycle is iterated until convergence.Let M (r) (β) denote the elastic net objective function with ℓ (r) (β) in place of ℓ(β).Let β(r,s) j denote the estimate of β j at the s th inner loop iteration of the r th outer loop iteration.Let M (r,s) j r) as a marginal function of the j th regression coefficient only, with all other coefficients fixed at their current estimates.The s th iteration coordinate descent update of the j th coefficient is arg min M (r,s) j (t).If c j = 0 (i.e.β j is unpenalized), then this can be solved straightforwardly by setting d dt M (r,s) j (t) = 0.In general, for t = 0, where X •j denotes the j th column of X, X •−j denotes X with the j th column deleted, and β(r,s) (t) runs from −∞ to +∞ and is monotonically increasing.The only point where the derivative does not exist is at t = 0, where it jumps up by 2λαc (t) = 0. Hence, the coordinate descent update can be written as β(r,s+1) where S(x, y) = sign(x)(|x| − y) + is the soft-thresholding operator, as defined in Friedman et al. (2010).On a side note, in some situations one may wish to include a model constraint that forces some or all of the coefficients to be nonnegative (e.g. to implement the nonnegative garrote penalty discussed in Breiman (1995)).With this constraint, the coordinate descent update becomes β(r,s+1) where the only difference is the restriction that t ≥ 0. When arg min M (r,s) The coordinate descent update can be written as

Computational efficiency and numerical stability
It would be computationally inefficient to compute the coordinate descent coefficient updates in the form presented in Section 3.5.First of all, W (r) is a sparse block diagonal matrix and should not be multiplied without taking this into consideration.Secondly, it is not necessary to explicitly compute z (r) , which depends on W −1 1 , W −1 2 , . . ., W −1 N .It is not necessarily computationally expensive to compute these K × K matrix inverses when K is small, but it is still better to avoid doing so.For computational efficiency, it is better to write the coordinate descent updates in terms of the score and information matrix.This is done by making the following two substitutions in (2) The subscripts on terms with square brackets indicate the j th vector element or matrix diagonal.Only the third term on the right hand side of the first substitution needs to be updated with each inner loop iteration because it depends on β(r,s) .The other terms are a function of β(r) , so they only need to be updated during the outer loop.
Furthermore, the X i matrices are sparse for ELMO models.In some programming languages, it may be advantageous to use this block structure to compute X i β and X T i W i X i .This has the additional advantage that it is not necessary to store X 1 , X 2 , . . ., X N in memory, but rather just x 1 , x 2 , . . ., x N , which are smaller.
Numerical instability of the information matrix can arise when the fitted class probabilities approach zero for any observation.This is problematic even if the near-zero probability occurs for an observation i and class j such that y ij = 0.A way to prevent numerical instability is to cap the fitted probabilities at some minimum value just for the information matrix calculation.(The pMin argument sets this threshold for the optimization function ordinalNet in the ordinalNet R package.)The score and likelihood are computed with uncapped fitted probabilities.

Regularization parameter sequence
Often we are interested in computing solutions for a sequence of λ values, rather than a single value.For α > 0, there always exists a threshold value λ max where the first coefficient enters the solution path.All penalized coefficients are set to zero for any λ > λ max .An off-the-shelf method to generate a reasonable sequence of λ values is to let λ min = 0.01 × λ max and consider a sequence of λ values that is uniform between λ max and λ min on the log scale (Friedman et al., 2010).
To calculate λ max , we first fit the intercept-only model by unpenalized maximum likelihood.Also include any unpenalized non-intercept coefficients if there are any.We then calculate the quadratic approximation at this solution.Each penalized coefficient has a threshold value of λ where its coordinate descent update becomes nonzero.The minimum threshold value among all coefficients is λ max , as this is the value where the first coefficient enters the solution path.Specifically, where β is the intercept-only maximum likelihood estimate, and W and z are calculated at β.

Starting values
Park and Hastie (2007) proposed an efficient estimation procedure for a decreasing sequence of λ values.The β solution for each λ value is used as the starting value for the next λ in the sequence.This technique is known as "warm starts."Furthermore, it is not necessary to update all coefficient estimates during the coordinate descent inner loop.Many coefficient estimates will begin and remain at zero throughout the inner loop, especially for larger λ values.It is more efficient to cycle through only the coefficients that had nonzero estimates at the previous step.The set of nonzero coefficients is known as the "active set."After the coordinate descent inner loop converges, one pass can be made over each coefficient outside the active set.If the coordinate descent update is zero for all of them, then the optimal solution has been reached.If the final pass changes any coefficient estimate from zero to a nonzero value, then the coordinate descent loop should be continued including these new nonzero coefficients in the active set.
A reasonable set of starting values can usually be obtained by passing the observed response category frequencies into the link function-this provides intercept starting values, and all other coefficients can start at zero.This is also the solution corresponding to λ max if there are no unpenalized non-intercept coefficients.If the first λ value is not λ max or some of the non-intercepts are unpenalized, then this is still usually a reasonable set of starting values for the first λ value.

Stopping rule
The coordinate descent procedure has an inner and outer loop, both of which require convergence definitions and thresholds.A suggestion is to define convergence using the relative change in the elastic net objective.For the outer loop, the definition is M( β(r) )−M( β(r−1) ) M( β(r−1) ) < ǫ out .For the inner loop, the quadratic approximation to the log-likelihood should be used instead of the true log-likelihood, so the definition is M (r) ( β(r,s) )−M (r) ( β(r,s−1) ) M (r) ( β(r,s−1) ) < ǫ in .A small constant can also be added to the denominator to allow convergence when the log-likelihood is near zero.Based on some trial and error, it seems efficient to set the outer and inner convergence thresholds, ǫ out and ǫ in , to the same value.

Issues with the cumulative probability family
The cumulative probability family has a constrained parameter space because the cumulative probabilities must be monotone for every x in the population.It was discussed in Section 2.5 that if the constraint is only enforced for x in the training sample, then difficulties may arise because an out-of-sample x may have fitted probabilities that are not monotone.
The parameter space constraint can also create difficulties during optimization.Although the likelihood is undefined outside the constrained parameter space, we could define an improper likelihood on the unrestricted parameter space.This improper likelihood would allow observations to have fitted probabilities greater than zero or less than one, and it would be defined as zero anywhere that an observation has negative probability in a class that was observed.This is essentially the likelihood that coordinate descent algorithm is designed to optimize.As a result, the algorithm may seek a solution that lies outside the constrained parameter space.
When the solution path leaves the constrained parameter space, we simply terminate the optimization algorithm at the λ value where this occurred.Further work is required to devise a constrained optimization procedure for the cumulative probability family.

Simulation
We have discussed three penalized ELMO model forms for ordinal data: the parallel model, the nonparallel model, and the semi-parallel model.The purpose of the following simulation experiments is to show scenarios where each of these three model types yields better out-of-sample prediction accuracy than the others.
We conducted three simulation experiments, each based on 100 replicates, i.e. simulated datasets.For each replicate, the data were simulated from a forward stopping ratio model with the parameters shown in Table 2.In all of the experiments, the covariates were simulated as independent, standard normal random variables.Now consider the estimation procedures.For each of the three models, parallel, nonparallel, and semi-parallel, the elastic net tuning parameter was set to α = 1 (i.e.lasso penalty).For the semi-parallel model, the tuning parameter ρ was set to one.A λ sequence of twenty values was generated uniformly on the log scale between λ max of the full training data and 0.01 × λ max .For each simulated dataset, five-fold cross validation was used to select the optimal tuning parameter.
Table 2: Training data parameters for the three simulation experiments.For each experiment, the data generating mechanism was a forward stopping ratio model.
Simulation results are shown in Table 3.For a given replicate, out-of-sample prediction accuracy was evaluated as the average log-likelihood of 10,000 observations generated from the same distribution as the training set.The validation data set was chosen large enough that the within-replicate standard error of the mean was negligible relative to the across-replicate standard error.More precisely, the results were produced in the following way.For a given replicate, generate 10,000 observations from the same distribution as the training set.Fix replicate i, for these 10,000 observations, compute log-likelihoods x i,j , j = 1, ..., 10000 and denote the sample average as v i , Table 3 reports the sample average of the v i 's (across the 100 simulated datasets), v = 1 100 100 i=1 v i and the standard error of v.For Simulation 1, the data generating model is nonparallel, and the sample size is relatively large.This results in the nonparallel and semi-parallel models having better fits than the parallel model.For Simulation 2, the data generating mechanism is parallel, and only one in three covariates has a nonzero effect.This results in the parallel and semi-parallel models having better fits than the nonparallel model.Simulation 3 is similar to Simulation 2, but the first covariate has a highly nonparallel effect.This results in the semi-parallel model having a better fit than both the parallel and nonparallel models.
Note that in every simulation scenario, the semi-parallel model fit is nearly as good, if not better, than both the parallel and nonparallel model fits.This is not to say that the semi-parallel model should always be preferred, but it is evidence that it is a highly versatile model.In practice, cross validation could be used to determine which of the three methods performs the best, using out-of-sample log-likelihood, or another measure of fit.In addition, it could be worthwhile to use cross validation to compare different values of ρ for the semi-parallel fit.

Method comparison
We now demonstrate ELMO class models alongside other methods for outof-sample prediction and variable selection.We use a dataset from a cancer genomics example presented by Archer et al. (2014).The data come from the Gene Expression Omnibus GSE18081.The CRAN package ordinalgmifs contains a filtered version of this dataset called hccmethyl.It contains expression levels of 46 genes from 56 human subjects.The measurements come from liver tissue samples assayed using the Illumina GoldenGat Methylation BeadArray Cancer Panel I (Archer et al., 2010).Twenty subjects have a normal liver, 16 have cirrhosis (disease), and 20 have hepatocellular carcinoma (severe disease).These categories have a natural ordering according to disease severity.
The analysis goal was to use gene expression values to predict liver diseasemore specifically, to estimate the probability that each subject's liver would be classified as healthy, diseased, or severely diseased.The number of predictors is large relative to the sample size, making regularization and/or variable selection imperative.We use five-fold cross validation to compare the out-of-sample prediction performance of various methods.We use out-of-sample log-likelihood and misclassification rate as performance criteria.A total of seven methods were compared.We summarize them below.
• λ was tuned by five-fold cross validation within each cross validation fold, selecting the value with the best average out-of-sample log-likelihood.
• λ was selected from the same sequence within each fold.λ max was calculated from the full data, and a sequence of twenty values was generated uniformly on the logarithmic scale from λ max to λ max /100.
• Two models: standard lasso and group lasso.
• λ was tuned by five-fold cross validation within each cross validation fold, selecting the value with the best average out-of-sample log-likelihood.
• λ was selected from the same sequence within each fold.λ max was calculated from the full data, and a sequence of twenty values was generated uniformly on the logarithmic scale from λ max to λ max /100.
• The solution path was fit with step size 0.01.
• AIC was used for tuning within each fold.Specifically, the model with the best AIC was selected from all of the models in the solution path.
4. Cumulative logit model with AIC, forward stepwise variable selection.
Results for the experiment are summarized in Table 4.The GMIFS method performed the best according to both performance criteria.The parallel cumulative logit lasso, the semi-parallel cumulative logit lasso, and the multinomial logistic regression with (standard) lasso performed comparably to the GMIFS method.The other methods did not perform well.
Poor performance of the nonparallel, cumulative logit lasso model can be attributed to the parameter space restriction discussed in Sections 2.5 and 3.11.With this dataset, it is easy for the nonparallel model to predict non-monotone cumulative probabilities for out-of-sample observations.However, the largest λ value in the sequence corresponds to the null model, which cannot have nonmonotone cumulative probabilities.As a result, the cross validation tuning procedure tends to select the null model or a highly penalized model; neither model is good for prediction.
The AIC stepwise method performed poorly in terms of out-of-sample loglikelihood.This poor performance is likely due to overfitting.More specifically, this method occasionally estimates a very small probability in the observed class of an out-of-sample observation.The out-of-sample misclassification rate is more reasonable than the out-of-sample log-likelihood, but other methods perform better.4: Out-of-sample log-likelihood and misclassification rates for the method comparison based on the liver disease dataset (GSE18081).The value reported is the mean (standard error) across five cross validation folds.

Log-likelihood
In addition to class prediction, penalized regression can be used for variable selection, i.e. to determine which of the 46 genes are most associated with liver disease.Table 5 shows which genes were selected by each method.All models were tuned on the full data the same way that they were trained in the cross validation study.
Table 5: Regression coefficient estimates for the method comparison based on the liver disease dataset (GSE18081).The column headers are abbreviations for the methods listed in the same order as Table 4.For the methods with multiple coefficients per predictor, the number of nonzero coefficient estimates is reported instead of the coefficient value.The nonparallel cumulative logit model can have up to two nonzero coefficients per predictor.The semi-parallel cumulative logit and multinomial logistic regression models can have up to three nonzero coefficients per predictor.By design, the group lasso penalty either selects all three coefficients or sets them all to zero.Note that the coefficient estimates from MASS::polr are multiplied by -1 to be consistent with the ordinalNet and ordinalgmifs parameterizations of the cumulative logit model.

The ordinalNet R package
The ordinalNet package contains the following functions: • ordinalNet is the main function for fitting parallel, nonparallel, and semiparallel regression models with the elastic net penalty.It returns an 'or-dinalNetFit' S3 object.
• ordinalNetTune uses K-fold cross validation to obtain out-of-sample loglikelihood and misclassification rates for a sequence of lambda values.
• ordinalNetCV uses K-fold cross validation to obtain out-of-sample loglikelihood and misclassification rates.Lambda is tuned within each cross validation fold.
Below is a description of the ordinalNet function and its arguments.

Arguments
x Covariate matrix.It is recommended that categorical covariates are converted to a set of indicator variables with a variable for each category (i.e.no baseline category); otherwise the choice of baseline category will affect the model fit.
y Response variable.Can be a factor, ordered factor, or a matrix where each row is a multinomial vector of counts.A weighted fit can be obtained using the matrix option, since the row sums are essentially observation weights.Non-integer matrix entries are allowed.
alpha The elastic net mixing parameter, with 0 <= alpha <= 1. alpha=1 corresponds to the lasso penalty, and alpha=0 corresponds to the ridge penalty.
standardize If standardize=TRUE, the predictor variables are scaled to have unit variance.Coefficient estimates are returned on the original scale.
penaltyFactors Nonnegative vector of penalty factors for each variable.This vector is multiplied by lambda to get the penalty for each variable.If NULL, the penalty factor is one for each coefficient.
positiveID Logical vector indicating whether each coefficient should be constrained to be non-negative.If NULL, the default value is FALSE for all coefficients.
family Specifies the type of model family.Options are "cumulative" for cumulative probability, "sratio" for stopping ratio, "cratio" for continuation ratio, and "acat" for adjacent category.
reverse Logical.If TRUE, then the "backward" form of the model is fit, i.e. the model is defined with response categories in reverse order.For example, the reverse cumulative model with K + 1 response categories applies the link function to the cumulative probabilities P (Y ≥ 2), . . ., P (Y ≥ K +1), rather then P (Y ≤ 1), . . ., P (Y ≤ K).
link Specifies the link function.The options supported are logit, probit, complementary log-log, and cauchit.Only used if customLink=NULL.
customLink Optional list containing a vectorized link function g, a vectorized inverse link h, and the Jacobian function of the inverse link getQ.The Jacobian should be defined as ∂h(η)/∂η T (as opposed to the transpose of this matrix).
parallelTerms Logical.If TRUE, then parallel coefficient terms will be included in the model.parallelTerms and nonparallelTerms cannot both be FALSE.
nonparallelTerms Logical.if TRUE, then nonparallel coefficient terms will be included in the model.parallelTerms and nonparallelTerms cannot both be FALSE.
parallelPenaltyFactor Numeric value greater than or equal to zero.Lambda is multiplied by this factor (as well as variable-specific penaltyFactors) to obtain the penalties for parallel terms.Only used if parallelTerms=TRUE.
lambdaVals An optional user-specified lambda sequence (vector).If NULL, a sequence will be generated based on nLambda and lambdaMinRatio.In this case, the maximum lambda is the smallest value that sets all penalized coefficients to zero, and the minimum lambda is the maximum value multiplied by the factor lambdaMinRatio.
nLambda Positive integer.The number of lambda values in the solution path.
Only used if lambdaVals=NULL.
lambdaMinRatio A factor greater than zero and less than one.Only used if lambdaVals=NULL.
includeLambda0 Logical.If TRUE, then zero is added to the end of the sequence of lambdaVals.This is not done by default because it can significantly increase computational time.An unpenalized saturated model may have infinite coefficient solutions, in which case the fitting algorithm will still terminate when the relative change in log-likelihood becomes small.Only used if lambdaVals=NULL.
alphaMin Value greater than zero, but much less than one.If alpha < alphaMin, then alphaMin is used to calculate the maximum lambda value.When alpha=0, the maximum lambda value would be infinite otherwise.
pMin Value greater than zero, but much less than one.During the optimization routine, the Fisher information is calculated using fitted probabilities.For this calculation, fitted probabilities are capped below by this value to prevent numerical instability.
stopThresh In the relative log-likelihood change between successive lambda values falls below this threshold, then the last model fit is used for all remaining lambda.
threshOut Convergence threshold for the coordinate descent outer loop.The optimization routine terminates when the relative change in the penalized log-likelihood between successive iterations falls below this threshold.It is recommended to set theshOut equal to threshIn.
threshIn Convergence threshold for the coordinate descent inner loop.Each iteration consists of a single loop through each coefficient.The inner loop terminates when the relative change in the penalized approximate loglikelihood between successive iterations falls below this threshold.It is recommended to set theshOut equal to threshIn.
maxiterOut Maximum number of outer loop iterations.
maxiterIn Maximum number of inner loop iterations.
printIter Logical.If TRUE, the optimization routine progress is printed to the terminal.
printBeta Logical.If TRUE, coefficient estimates are printed after each coordinate descent outer loop iteration.
warn Logical.If TRUE, the following warning message is displayed when fitting a cumulative probability model with nonparallelTerms=TRUE (i.e.nonparallel or semi-parallel model)."Warning message: For out-of-sample data, the cumulative probability model with nonparallelTerms=TRUE may predict cumulative probabilities that are not monotone increasing." The warning is displayed by default, but the user may wish to disable it.
keepTrainingData Logical.If TRUE, then x and y are saved with the returned "ordinalNetFit" object.This allows predict.ordinalNetFit to return fitted values for the training data without passing a newx argument.

Demonstration in R
This section contains six examples that demonstrate different aspects of the ordinalNet package.Specifically, using the GSE18081 dataset from the Gene Expression Omnibus, we demonstrate how the penalized ELMO class models can be used for prediction and variable selection.

Example 1
We illustrate how to fit the Gene Expression Omnibus GSE18081 data using ordinalNet.We fit a cumulative probability model with logit link (proportional odds model).We use the default settings parallelTerms=TRUE and nonparallelTerms=FALSE to fit the parallel model.We use the default elastic net tuning parameter alpha=1 to select the lasso penalty.We use the default settings of lambdaVals=NULL, nLambda=20, and lambdaMinRatio=1e-2 to generate a sequence of twenty λ values, with λ max equal to the smallest value that sets every coefficient to zero and λ min = λ max • 0.01.The sequence runs from λ min to λ max uniformly on the logarithmic scale.The summary method displays the lambda sequence (lambdaVals), number of nonzero coefficients (nNonzero), the log-likelihood (loglik), percent deviance explained (pctDev), and AIC and BIC.The AIC and BIC are calculated using nNonzero as the approximate degrees of freedom.The coef method returns the coefficient estimates of any model fit in the sequence-the best AIC fit is selected by default.The matrix=TRUE option returns the coefficients in matrix form with a column corresponding to each linear predictor.Because fit1 is a parallel model, the coefficient columns are identical except for the intercepts.

Example 5
The ordinalNetCV function uses K-fold cross validation to obtain out-of-sample performance of models that are tuned within each cross validation fold.We use the default value of nFolds=5 and the same default λ sequence as above.We also use the default settings of nFoldsCV=5 and tuneMethod="cvLoglik" to tune λ by 5-fold cross validation within each fold, each time selecting the λ value with the best average out-of-sample log-likelihood.
We compare the performance of the parallel, semi-parallel, and nonparallel cumulative probability models.By default, the ordinalNetCV function will randomly create the fold partitions, but we pass a list of fold partitions to ensure that the same partitions are used for each of the three methods.(This could also be done by setting the same seed before each of the three calls to ordinalNetCV).Although not of critical importance, we also increase the outer iteration limit to 200 for the semi-parallel model.
We see that the parallel and semi-parallel models have similar performance, but the nonparallel model is much worse.This is because the nonparallel model out-of-sample log-likelihood becomes undefined (non-monotone cumulative probabilities) within the first few λ values.Further examination of the λ index value selected for each fold reveals that often the first λ value (i.e. the null model) is selected for the nonparallel model.It likely does not matter that the solution path is terminated for leaving the parameter space because the cross validation procedure would probably select very large λ values even if the entire solution path were available.

Discussion
This paper introduced the elmentwise link, multinomial-ordinal (ELMO) model class, a rich class of multinomial regression models that includes commonly used categorical regression models.Each of these models has both a parallel and nonparallel form.The parallel form is appropriate for ordinal data, while the nonparallel form is a more flexible model which can be used with an unordered categorical response.We also introduced the semi-parallel model form, which can be used with the elastic net penalty to shrink the nonparallel model toward the parallel model.
The motivation for this work began with a need to extend variable selection tools for ordinal logistic regression.For instance, consider the problem of developing a gene signature to predict response to a novel therapy, where the observed patient response belongs to one of the following categories: no response, partial response, or complete response.We developed these tools in the general ELMO framework.Specifically, we proposed a coordinate descent fitting algorithm for the ELMO class with the elastic net penalty.The algorithm is general and can also be applied to multinomial regression models outside the ELMO class.
We considered numerical experiments to highlight different features of the model class and to demonstrate the use of the related R code.We presented different simulation scenarios to demonstrate cases where the parallel, nonparallel, and semi-parallel each achieved better out-of-sample prediction performance than the other two models.With the Gene Expression Omnibus GSE18081 data set, we demonstrated the use of the penalized ELMO class for prediction and variable selection.Finally, we introduced the R package ordinalNet, which implements a coordinate descent algorithm for parallel, nonparallel, and semi-parallel models of the ELMO class.It is available on the Comprehensive R Archive Network at http://CRAN.R-project.org/package=ordinalNet.
We consider two possible directions for future research: code speedup via C ++ and questions of statistical inference.Rcpp and RcppArmadillo are R packages which allow integration of C ++ code into R (Eddelbuettel et al., 2011;Eddelbuettel, 2013;Eddelbuettel and Sanderson, 2014).Our code is written with separate functions for the inner and outer coordinate descent loops.Because of the number of calls to it in a typical run of the algorithm, the inner loop, in particular, could benefit from speed up via C ++ .
The ordinalNet package does not provide standard errors for estimates.We quote a relevant section from the penalized vignette (Goeman et al., 2014).
It is a very natural question to ask for standard errors of regression coefficients or other estimated quantities.In principle such standard errors can easily be calculated, e.g. using the bootstrap.Still, this package deliberately does not provide them.The reason for this is that standard errors are not very meaningful for strongly biased estimates such as arise from penalized estimation methods.Penalized estimation is a procedure that reduces the variance of estimators by introducing substantial bias.The bias of each estimator is therefore a major component of its mean squared error, whereas its variance may contribute only a small part.
The topic of post-selection inference has been studied in both the classic setting (Zhang, 1992;Leeb and Pötscher, 2005;Wang and Lagakos, 2009;Berk et al., 2013), where the number of observations exceeds the number of predictors, and the high-dimensional setting (Javanmard and Montanari, 2014;Lockhart et al., 2014;Tibshirani et al., 2016a).In the high-dimensional setting, we would like to highlight the groundbreaking work of Lockhart et al. (2014).They proved the asymptotic distribution of their test statistic specifically for the linear model, but their simulation results suggest that the same test statistic could be used for generalized linear models.This work may provide a path for post-selection inference for penalized multinomial and ordinal regression models.A Uniqueness of the semi-parallel model estimator This section addresses the uniqueness problem for the semi-parallel model.The semi-parallel model without the elastic net penalty is not identifiable, as there are infinitely many parametrizations for any particular model.However, different parametrizations have different elastic net penalty terms; therefore, the penalized likelihood favors some parametrizations over others.We will demonstrate that amongst almost all parametrizations for a given model, the elastic net penalty has a unique optimum; hence, the penalized likelihood has a unique optimum.There is one exception: the lasso penalty with integer-valued ρ.In this case, there may be a small range of optima on a closed interval.We proceed in the following manner.First, we formulate the basic problem.Next we consider uniqueness for the ridge penalty (α = 0), which is the simplest case.We then consider the lasso penalty (α = 1), where this exception can occur.Finally, we consider the elastic net penalty with α ∈ (0, 1).These derivations are related to derivations for the lasso in the linear model setting (Osborne et al., 2000;Rinaldo, 2008;Tibshirani, 2013).
To see this, for any ζ set δ k = β k − ζ for all k.All of these parametrizations have the same likelihood because they specify the same model, but they have different elastic net penalty terms proportional to Our goal is to find the value of ζ that minimizes the elastic net penalty.

Ridge regression
In this case, the Lagrangian is differentiable everywhere.Consider Solving this yields Solving this yields The solution is unique for any ρ ≥ 0. Next, consider

Lasso
We want the solution where ∂L ∂ζ equals or crosses zero.That is, we are searching for the value of ζ where f (ζ) equals or crosses ρ, with Note that if ρ ≥ K, then the solution will be ζ = 0. Hence, if ρ ≥ K, then all parallel coefficients will be penalized to zero, and the fit will be equivalent to the nonparallel model with the elastic net penalty.Now, if ρ is not an integer, then the solution is unique.If ρ is an integer, then the solution could be unique, or there may be a range of solutions on a closed interval between two consecutive β's, ranked by value.Solving this yields We want the solution where ∂L ∂ζ equals or crosses zero.This solution is less transparent than that of ridge or lasso.However, the partial derivative is piecewise linear in ζ and never constant over a range of values.Hence, the solution is unique.

B The inverse link function and its Jacobian for specific MO families
Each MO family is defined by the g MO (p) component of its multivariate link function.For optimization, it is necessary to compute the inverse h MO (δ) and its Jacobian Dh MO (δ).In this section, we provide a method to compute these three functions for the cumulative probability, stopping ratio, continuation ratio, and adjacent category families.In some cases, it is convenient to compute the elements of these functions recursively.Although the Jacobian is, strictly speaking, a function of δ, we write it in terms of both p and δ when convenient.This can be done because there is a one-to-one correspondence between δ and p.
Each family has a forward and backward form.We present only one of these forms for each family.To fit the backward form, one can simply define the response categories in reverse order and fit the forward model, and vice versa.

Forward adjacent category family
This family is defined by δ j = P (Y = j + 1|j ≤ Y ≤ j + 1) (see Table 1).From this definition, for all j, we have [g MO (p)] j = p j+1 p j + p j+1 .

C Quadratic approximation to the log-likelihood
The quadratic approximation ℓ (r) (β) is the second order Taylor expansion of ℓ at β(r) , up to an additive constant that does not depend on β.Also, in ℓ (r) (β), the Hessian is replaced by its expectation, −I( β(r) ).We give the derivation below, where  Goeman et al. (2014) assert that while standard errors of penalized regression coefficients could be calculated, e.g. by the bootstrap, these standard errors should be interpreted with care.In particular, they should not be used to construct confidence intervals.This is because penalized regression coefficient estimators have substantial bias.For this reason, they deliberately do not provide standard error estimates in the penalized package.In the classical setting, where the number of observations exceeds the number of predictors, then one can use the usual likelihood-based methods to construct confidence intervals and hypothesis tests.However, even in this setting, such methods can lead to incorrect inference (Zhang, 1992;Leeb and Pötscher, 2005;Wang and Lagakos, 2009;Berk et al., 2013).Tibshirani and Taylor (2011), Lockhart et al. (2014), andTibshirani et al. (2016b) laid the groundwork of an asymptotic inference method for penalized regression.They proved the asymptotic distribution of their test statistic specifically for least squares regression, but their simulation results suggest that the same test statistic could be used for generalized linear models.This may be a path toward high dimensional inference for penalized multinomial and ordinal regression models.

Table 1 :
Examples of multinomial-ordinal (MO) families.For each example, Y is a categorical random variable with class probability vector (p 1 , p 2 , . . ., p K then the derivative attains zero, and the value at which this occurs is arg min M By setting parallelTerms=TRUE and nonparallelTerms=TRUE, we obtain the semi-parallel model fit for fit2.Because this is a cumulative probability model, we set warn=FALSE to suppress the warning that the semi-parallel form is susceptible to non-monotone cumulative probabilities for out-of-sample predictions.We use the default semi-parallel tuning parameter ρ, which is parallelPenaltyFactor=1.The coefficient matrix of the best AIC fit has nearly identical columns for the first several predictors, but they differ for the first predictor.We now demonstrate the problem with the nonparallel cumulative probability model discussed in Sections 2.5 and 3.11.(Thesemi-parallel model is also susceptible to this issue, but it is less prone.It does not occur for the semiparallel model on this data set).As seen from the summary method output, the solution path is terminated after the third λ value where the optimum leaves the constrained parameter space.Better solutions must exist for the remaining λ values, but the coordinate descent optimization procedure is not designed handle this issue.The ordinalNetTune function uses K-fold cross validation to obtain out-ofsample performance for a sequence on λ values.We use the default setting of nFolds=5 and the default sequence of twenty λ values obtained from the model fit to the full data.The user can use this information to tune the model, for example by selecting the λ value with the best average out-of-sample likelihood across folds, as demonstrated below.