ordinalgmifs: An R Package for Ordinal Regression in High-dimensional Data Settings

Abstract

High-throughput genomic assays are performed using tissue samples with the goal of classifying the samples as normal < pre-malignant < malignant or by stage of cancer using a small set of molecular features. In such cases, molecular features monotonically associated with the ordinal response may be important to disease development; that is, an increase in the phenotypic level (stage of cancer) may be mechanistically linked through a monotonic association with gene expression or methylation levels. Though traditional ordinal response modeling methods exist, they assume independence among the predictor variables and require the number of samples (n) to exceed the number of covariates (P) included in the model. In this paper, we describe our ordinalgmifs R package, available from the Comprehensive R Archive Network, which can fit a variety of ordinal response models when the number of predictors (P) exceeds the sample size (n). R code illustrating usage is also provided.

Keywords

ordinal response high-dimensional features penalized models R

Introduction

Though there are traditional methods for modeling an ordinal response, these methods assume independence among the predictor variables and require that the number of samples (n) exceed the number of covariates (P) included in the model. For high-throughput genomic data, P >> n, so the traditional ordinal response models cannot be estimated. Penalized models have been demonstrated to have excellent performance when applied to high-throughput genomic datasets in fitting linear, logistic, and Cox proportional hazards models. However, penalized methods have not been fully extended to the ordinal response setting even though ordinal responses are pervasive in medicine. For example, tissue samples may be collected with the goal of classifying them as normal < pre-malignant < malignant or by stage of cancer. In these cases, molecular features monotonically associated with the ordinal response may be important to disease development; that is, an increase in the phenotypic level (stage of cancer) may be mechanistically linked through a monotonic association with gene expression. While one can apply nominal response classification methods to ordinal response data, in so doing some information is lost that may improve the predictive performance of the classifier. We previously developed software for the R programming environment that uses the coordinate descent algorithms of Refs. 1 and 2 for fitting penalized constrained continuation ratio models.³ These packages, glmpathcr and glmnetcr, are capable of modeling an ordinal response in settings where P >> n. Unfortunately, extension to the cumulative link, adjacent category, stereotype logit models, and other link functions using this framework is not straightforward. Other algorithms can be used for obtaining solutions for the least absolute shrinkage and selection operator (LASSO)^4,5 and elastic net penalized models.⁶ In the linear regression setting, the incremental forward stagewise (IFS) is a penalized solution that enforces monotonicity.⁷ IFS can be generalized to problems involving those other than squared error loss, and the adaption is called the generalized monotone incremental forward stagewise (GMIFS) method.⁷ Herein, we extended the GMIFS method⁷ to ordinal response models and describe our ordinalgmifs R package. The ordinal.gmifs function can be used to fit traditional and penalized cumulative link, forward continuation ratio, and backward continuation ratio models using either a logit, probit, or complementary log–log link. It can also be used to fit adjacent category and stereotype logit models.

In the following sections, we provide an overview of the LASSO, IFS, and GMIFS methods for the linear and logistic regression settings. We then describe our implementation of the GMIFS method for modeling an ordinal response. Our implementation includes methods for initializing the ordinal thresholds, methods for updating estimates for a no penalty subset of predictors if specified by the user, derivatives for identifying which covariate to update, and convergence criteria. Because the GMIFS method requires the derivatives of each log-likelihood, where the log-likelihood function is derived from the traditional ordinal response models, including the cumulative logit, adjacent category, continuation ratio, and stereotype logit models, we also provide an overview of traditional ordinal response models. We then illustrate the functions in the ordinalgmifs R package using a dataset where we were interested in predicting normal < pre-neoplastic < neoplastic states of liver disease using high-throughput methylation data.⁸

Brief Overview of Penalized Methods

Traditional variable selection methodologies, including best subsets, and forward selection and backward elimination procedures, often fail to provide feasible and stable results because of their strong assumptions of covariate independence as well as their discrete procedures. Regularization methods provide a more continuous process and yield statistical models having coefficients with non-zero estimates for “important” covariates, while many coefficients are shrunk to be exactly zero. There are several algorithms that yield a penalized solution that are applicable for continuous responses.

Penalized Methods for Linear Regression.

The LASSO⁴ is a widely used method for deriving a parsimonious model for high-dimensional data. It produces penalized estimates of the unknown regression parameters by including the L₁ norm of the coefficients as a constraint in the least-squares estimate

{\hat{β}}^{lasso} = \underset{β}{\arg \min} (\sum_{i = 1}^{n} {(y_{i} - β_{0} - \sum_{p = 1}^{P} x_{i j} β_{p})}^{2} + λ \sum_{p = 1}^{P} | β_{p} |);

hence, the terms LASSO and L₁ penalized regression are often used synonomously. Here λ is the tuning parameter that controls the amount of shrinkage: if λ = 0, the solution is the OLS estimates; as λ increases, the amount of shrinkage of the parameter estimates increases. Because the L₁ model shrinks some coefficients to be exactly zero, it is better than ridge regression in terms of model parsimony and interpretability. Algorithms for fitting penalized linear regression models include the least angle regression (LAR) algorithm⁹ and the IFS regression method.⁷ The advantage of IFS is that it can be generalized to problems involving other than squared error loss. Therefore, we present the IFS method for the linear regression setting as a prelude to the GMIFS method.

IFS Method for Linear Regression.

The forward stagewise method is a greedy procedure similar to forward stepwise regression but is more cautious in that the coefficient updates are made using very small increments. It was historically considered to be inefficient because of its slow-fitting nature but later showed to be competitive in terms of variable selection stability and prediction accuracy.⁷ The IFS method⁷ proceeds according to the following steps to admit the solution: 1.

Standardize the predictors, and then initialize the residuals to $r = y - \bar{y}$ and initialize the p = 1, …, P parameter estimates to ${\hat{β}}_{1}, …, {\hat{β}}_{P} = 0$ .

Find the predictor x_m most correlated with r, $m = \underset{p}{\arg \max} (| \hat{ρ} (r, x_{P}) |)$ .

Update ${\hat{β}}_{m} = {\hat{β}}_{m} + δ_{m}$ , where $δ_{m} = ε \times s i g n (\hat{ρ} (r, x_{m}))$ and ε is some small positive constant such as 0.001.

Update r = r − δ_mx_m, and repeat steps 2–3 until no predictor has any correlation with r.

At the end of this iterative algorithm, the final ${\hat{β}}_{1}, …, {\hat{β}}_{P}$ are taken as the penalized solution. It was previously shown that the penalized solution can be obtained using the expanded predictor space $\hat{X} = [X : - X]$ .⁷ When using the expanded predictor space, the absolute value of the correlation estimates is unnecessary in step 2, while the coefficient and residual updates in steps 3 and 4 let δ_m = ε. Finally, the P × 1 vector of coefficient estimates corresponding to the original representation of the covariate matrix X can be obtained using ${\hat{β}}_{p} = {\hat{β}}_{p} - {\hat{β}}_{p + P}$ for p = 1, …, P.

Extending IFs for Logistic Regression.

For observations i = 1, …, n, let y_i represent the dichotomous dependent variable and

π (x_{i}) = P (y_{i} = 1 | x_{i}) = \frac{\exp (α + x_{i}^{T} β)}{1 + \exp (α + x_{i}^{T} β)}

represent the conditional probability of experiencing the event given the independent predictor variables x_i. Then the likelihood function is simply the product of the n independent binomial terms

L (α, β | y, x) = \prod_{i = 1}^{n} π {(x_{i})}^{y_{i}} {(1 - π (x_{i}))}^{1 - y_{i}} .

Mathematically, it is easier to maximize the log-likelihood, which is given by

\log (L (α, β | y, x)) = \sum_{i = 1}^{n} (y_{i} \log (π (x_{i})) + (1 - y_{i}) \log (1 - π (x_{i}))) .

When the response is discrete, minimizing the residual sum of squares is not reasonable; so a modified penalized model that maximizes the penalized log-likelihood is needed. For an L₁ penalized logistic regression model, this is expressed as

\log (L (α, β | y, x)) - {λ \sum}_{p = 1}^{P} | β_{p} | .

Obviously, the residual, r, in the IFS procedure is no longer appropriate for a categorical response, so it cannot be used in the logistic regression setting. However, the GMIFS regression method, where the coefficient to be updated is selected based on the gradient of the log-likelihood, can be used to provide the coefficient estimates for a penalized logistic regression model.⁷ Starting with the expanded covariate space $\tilde{X}$ , the GMIFS algorithm is 1.

Standardize the predictors and initialize the p = 1, …, 2P parameter estimates to ${\hat{β}}_{1}, …, {\hat{β}}_{P}, {\hat{β}}_{p + 1}, …, {\hat{β}}_{2 P} = 0$ at step s = 0.

Find the predictor x_m, where $m = \underset{p}{\arg \min} - δ \log L / δ β_{p}$ at the current estimate $\hat{β} = {\hat{β}}^{(s)}$ .

Update the corresponding coefficient ${\hat{β}}_{m}^{s + 1} = {\hat{β}}_{m}^{(s)} + ε$ to yield a new vector of parameter estimates ${\hat{β}}^{s + 1}$ .

Repeat steps 2 and 3 many times.

At the end of the iterative procedure, the P × 1 vector of coefficient estimates is obtained by ${\hat{β}}_{p} = {\hat{β}}_{p} - {\hat{β}}_{p + P}$ for p = 1, …, P.

Proposed GMIFS Algorithm for Ordinal Response Models.

Although the GMIFS algorithm has been described, there was no specific stopping criteria, but rather step 4 in Ref. 7 recommends to repeat steps 2 and 3 many times. Therefore, we implemented the criteria to stop the iterative process when the difference between two successive log-likelihoods is smaller than a pre-specified tolerance τ. Further, recall that for the linear regression scenario, the constant term β₀ term is commonly omitted because the response is centered. However, for ordinal response models having K different ordinal levels, we require K − 1 separate threshold estimates be included in the model. These K − 1 α intercept terms are essential for an ordinal response model that assumes proportional odds, because under this assumption, the β estimates do not have category-specific effects.¹⁰ Therefore, the log-odds ratio, which measures the degree of association between two ordinal levels, can be explained by the difference between the corresponding thresholds, logit(γ_j) − logit(γ_j’) = α_j − α_j’, such that the α estimates are highly important for distinguishing between ordinal classes. Moreover, there are some research problems for which a subset of predictor variables is to be included in the model (ie, not penalized). That is, in certain modeling contexts, it is desired to include relevant variables based on subject-specific knowledge. Therefore, we have included in our algorithm the ability to include a set of covariates as a no penalty subset of predictor variables. For K ordinal classes and P predictors, our GMIFS algorithm is: 0a.

Let $\tilde{X} = [X : - X]$ be our standardized predictors that are to be penalized and, if applicable, W be our predictor(s) that is(are) not to be penalized (no penalty subset). Let β represent the coefficients associated with X and θ represent the coefficients associated with W.

0b.

Initialize α and θ based on the specific ordinal response model (see the following ordinal response models section).

Initialize the components of ${\hat{β}}^{(s)}$ at step s = 0 as ${\hat{β}}_{1} = {\hat{β}}_{2} = … = {\hat{β}}_{P} = {\hat{β}}_{P + 1} = … = {\hat{β}}_{2 P} = 0$ .

Find $m = \underset{p}{\arg \min} - δ \log L / δ β_{p}$ at the current estimate ${\hat{β}}^{(s)}$ .

Update ${\hat{β}}_{m}^{(s + 1)} \leftarrow {\hat{β}}_{m}^{(s)} + ε$ .

Update α's and θ by maximum likelihood considering the vector of coefficients ${\hat{β}}^{s + 1}$ from step 3 as fixed.

Repeat steps 2–4 until log L^(s+1) − log L^(s) < τ, where τ is pre-specified tolerance.

Therefore, for each ordinal response model, we need to initialize the α₁, …, α_K-1 and θ estimates, provide expressions for δ log L/δβ_p for updating the β's, and provide maximum likelihood equations for updating the α1, …, α_K-1 and θ estimates after the β's have been updated. These equations are presented in the subsequent sections for the cumulative link, adjacent category, forward and backward continuation ratio model, and stereotype logit model. While our ordinalgmifs R package additionally allows fitting cumulative link, and forward and backward continuation ratio models using a probit or complementary log-log link, we have omitted the details regarding their derivatives and initial α estimates here.

Ordinal Response Models

Let Y_i represent the ordinal response for observation i that can take on one of K ordinal levels. Denote the n × P covariate matrix as x so that x_i represents a P × 1 vector for observation i and x_p represents the n × 1 vector for covariate p. For observations i = 1, …, n, the response Y_i can be reformatted as a response matrix consisting of n rows and K columns where

y_{i j} = {\begin{matrix} \begin{matrix} 1 & i f o b s e r v a t i o n i i s c l a s s j, \\ 0 & o t h e r w i s e . \end{matrix} \end{matrix}

Therefore, y_j is an n × 1 vector representing class j membership. Letting π_j(x_i) represent the probability that observation i with covariates x_i belongs to class j, the likelihood for an ordinal response model with K ordinal levels can be expressed as

L = \prod_{i = 1}^{n} \prod_{j - 1}^{K} π_{j} {(x_{i})}^{y_{i j}} .

(1)

Similar to logistic regression, it is more convenient to work with the log-likelihood, which be expressed as

\log L = \sum_{i = 1}^{n} \sum_{j = 1}^{K} y_{i j} \log (π_{j} (x_{i})) .

(2)

Cumulative Logit Model.

The cumulative logit model models K − 1 logits of the form.

P (Y_{i} \leq j | x_{i}) = \frac{\exp (α_{j} + x_{i}^{T} β)}{1 + \exp (α_{j} + x_{i}^{T} β)},

where α_j denotes the class-specific intercept and β is a p × 1 vector of coefficients associated with explanatory variables x_i.¹¹ Note that the class-specific probabilities can be calculated by subtracting successive cumulative logits,

π_{j} (x_{i}) = P (Y_{i} \leq j | x_{i}) - P (Y_{i} \leq j - 1 | x_{i}) .

Therefore, for any class j, we can express the class-specific probabilities by

π_{j} (x_{i}) = \frac{\exp (α_{j} + x_{i}^{T} β)}{1 + \exp (α_{j} + x_{i}^{T} β)} - \frac{\exp (α_{j - 1} + x_{i}^{T} β)}{1 + \exp (α_{j - 1} + x_{i}^{T} β)} .

Because we require α0 = − ∞, for j = 1, this expression simplifies to

π_{1} (x_{i}) = \frac{\exp (α_{1} + x_{i}^{T} β)}{1 + \exp (α_{1} + x_{i}^{T} β)},

and because we require α_K = ∞, for j = K, this expression simplifies to

π_{K} (x_{i}) = 1 - \frac{\exp (α_{K - 1} + x_{i}^{T} β)}{1 + \exp (α_{K - 1} + x_{i}^{T} β)} .

Therefore, the derivative of the log-likelihood with respect to βp is

\begin{matrix} \frac{δ \log L}{δ β_{p}} = x_{p}^{T} (\frac{y_{1}}{1 + \exp (α_{1} + x β)} \\ - \sum_{j = 2}^{K - 1} \frac{(\exp (α_{j} + α_{j - 1} + 2 x β) - 1) y_{j}}{(1 + \exp (α_{j} + x β)) (1 + \exp (α_{j - 1} + x β))} \\ - \frac{\exp (α_{K - 1} + x β) y_{K}}{1 + \exp (α_{K - 1} + x β)}) . \end{matrix}

(3)

In the GMIFS algorithm, the alpha terms are initialized to $α_{j} = \log i t (\sum_{i = 1}^{n} \sum_{k = 1}^{j} y_{i k} / n)$ . If a no penalty subset is included, we can consider x to be the expanded matrix of w: x, and the parameter vector β to be the expanded vector θ: β for calculating the class-specific probabilities. However, the derivative is only estimated with respect to the variables in the penalized portion of the model. Also, when a no penalty subset is included, α and θ are estimated by maximum likelihood at the initial step and after each β update. Because the class-specific probabilities are obtained by subtracting successive cumulative logits, we require α₁ ≤ α₂ ≤ ··· ≤ α_K-1. Therefore, the maximum likelihood estimation of θ and θ includes that constraint on the threshold estimates.

Adjacent Category Model.

The adjacent category model models the logits of ordinal categories that are adjacent to one another. That is, for any ordinal level j = 1, …, K − 1, its logit is

\log i t (γ_{j}) = α_{j} + x β,

(4)

where

γ_{i j} = \frac{π_{i, j + 1}}{π_{i j} + π_{i, j + 1}} .

(5)

Plugging equation 5 into the logit in equation 4, we have the simplification,

\begin{matrix} \log i t (γ_{i j}) = \log (\frac{\frac{π_{i, j + 1}}{π_{i j} + π_{i, j + 1}}}{1 - \frac{π_{i, j + 1}}{π_{i j} + π_{i, j + 1}}}) \\ = \log (\frac{π_{i, j + 1}}{π_{i j}}) \end{matrix}

so that

\log (\frac{π_{i, j + 1}}{π_{i j}}) = α_{j} + x_{i}^{T} β .

The class-specific probabilities for the adjacent category model can be calculated using the baseline category framework. Suppose we let the first class be our baseline category (note that the baseline category is arbitrary).

Then for any other class j ≤ K, we can express its baseline category logit as

\begin{matrix} \log (\frac{π_{j}}{π_{1}}) = \log (\frac{π_{j}}{π_{j - 1}}) + \log (\frac{π_{j - 1}}{π_{j - 2}}) + … + \log (\frac{π_{2}}{π_{1}}) \\ = α_{j - 1} + x β + α_{j - 2} + x β + … + α_{1} + x β \\ = \sum_{k = 1}^{j - 1} α_{k} + (j - 1) x β . \end{matrix}

Adding log(π₁) to both sides of this equation and then exponentiating yields

π_{j} = \exp (\sum_{k = 1}^{j - 1} α_{k} + (j - 1) x β + \log (π_{1}))

and

π_{1} = \frac{1}{1 + \sum_{j = 1}^{K - 1} \exp (\sum_{k = 1}^{j} α_{k} + j x β)} .

The derivative of the log-likelihood with respect to βp is the sum of the derivatives of the K class components

\frac{δ \log L}{δ β_{p}} = - x_{p}^{T} (\sum_{j = 1}^{K} [(j - 1) y_{j} - y_{j} (\frac{\sum_{l = 1}^{K - 1} l \exp (\sum_{k = 1}^{l} α_{k} + k x β)}{1 + \sum_{l = 1}^{K - 1} \exp (\sum_{k = 1}^{l} α_{k} + k x β)})) .

In the GMIFS algorithm, the α terms are initialized as the logits of the adjacent class probabilities: α₁ = log(π₂/π₁), α₂ = log (π₃/π₂), α_K-1 = log(π_K/π_K-1). If a no penalty subset is included, α and θ are estimated by maximum likelihood at the initial step and after each β update.

Continuation Ratio Models.

Although the likelihood for continuation ratio models can be expressed using equation 1, it is more convenient to express the likelihood using the continuation ratios. There are different ways in which one can set up a continuation ratio model, so here we present the backward and forward continuation ratio models separately.

Backward Continuation Ratio Model.

The backward continuation ratio model models the logit of the j = 2, …, K conditional probabilities or

\log i t (P (Y = j | Y \leq j, x)) = \log (\frac{P (Y = j | Y \leq j, x)}{P (Y < j | Y \leq j, x)}) = α_{j} + x β .

Here we have used the backward formulation, which is commonly used when progression through disease states from none, mild, moderate, and severe is represented by increasing integer values, and interest lies in estimating the odds of more severe disease compared to less severe disease.¹² Letting δ_ij represent the conditional probabilities,

δ_{i j} = δ_{j} (x_{i}) = P (Y_{i} = j | Y_{i} \leq j, x_{i}) = \frac{\exp (α_{j} + x_{i}^{T} β)}{1 + \exp (α_{j} + x_{i}^{T} β)},

such that for K ordinal classes, there are K − 1 logits. The likelihood can be expressed using these j = 2, …, K conditionally independent probabilities:

L (β | y, x) = \prod_{i = 1}^{n} \prod_{j = 2}^{K} δ_{i j}^{y_{i j}} {(1 - δ_{i j})}^{1 - \sum_{k = j}^{K} y_{i k}},

which can be factored as the product of K − 1 binomial likelihoods. Using this expression, the derivative of the log-likelihood with respect to β_p for the backward continuation ratio is given by

\frac{δ \log L}{δ β_{p}} = \sum_{j = 2}^{K} x_{p}^{T} (\frac{y_{j}}{1 + \exp (α_{j} + x β)} + \frac{(\sum_{k = j}^{K} y_{k} - 1) \exp (α_{j} + x β)}{1 + \exp (α_{j} + x β)}) .

(7)

In the GMIFS algorithm, the α terms are initialized as α_j = log (π/(1 − π_j)), where $Π_{j} = \sum_{i = 1}^{n} y_{i j} / \sum_{i = 1}^{n} \sum_{k = 1}^{j} y_{i k} .$ . If a no penalty subset is included, α and θ are estimated by maximum likelihood at the initial step and after each β update.

Forward Continuation Ratio Model.

The forward continuation ratio model models the logit of the conditional j = 1, …, K − 1 probabilities or

\log i t (P (Y = j | Y \geq j, x)) = \log (\frac{P (Y = j | Y \geq j, x)}{P (Y > j | Y \geq j, x)}) = α_{j} + x β .

As with the backward continuation ratio model, the likelihood for the forward continuation ratio model can be expressed using the K − 1 conditionally independent probabilities,

L (β | y, x) = \prod_{i = 1}^{n} {\prod_{j = 1}^{K - 1} δ_{i j}^{y_{i j}} (1 - δ_{i j})}^{\sum_{k = j}^{K} y_{i k} - y_{i j}} .

The derivative of the log-likelihood is the same, which for the forward continuation ratio model is given by,

\frac{δ \log L}{δ β_{p}} = \sum_{j = 1}^{K - 1} x_{p}^{T} (\frac{y_{j}}{1 + \exp (α_{j} + x β)} - \frac{(\sum_{k = j}^{K} y_{k} - y_{j}) \exp (α_{j} + x β)}{1 + \exp (α_{j} + x β)}) .

(8)

In the GMIFS algorithm the α terms are initialized as α_j = log (π_j/(1 − π_j)), where $Π_{j} = \sum_{i = 1}^{n} y_{i j} / \sum_{i = 1}^{n} \sum_{k = j}^{K} y_{i k} .$ . If a no penalty subset is included, α and θ are estimated by maximum likelihood at the initial step and after each β update.

Stereotype Logit Model.

The proportional odds versions of the cumulative logit, adjacent category, and forward and backward continuation ratio models all assume that the odds at x = x₁ compared to x = x₂ are exp(β(x₁ − x₂)) regardless of the response category. The stereotype logit model is more flexible than the proportional odds versions of the cumulative logit, adjacent category, forward continuation ratio, and backward continuation ratio models, yet is more parsimonious in comparison to fully unconstrained versions of these models.^13,14 From the definition of the stereotype logit model,¹³ the probability that observation i is from class j is defined as

π_{j} (x_{i}) = \frac{\exp (α_{j} + ϕ_{j} x_{i}^{T} β)}{\sum_{h = 1}^{K} \exp (α_{h} + ϕ_{h} x_{i}^{T} β)},

(9)

and to impose an ordinality restriction, we let α_K = 0 and ϕ_K = 0 such that for category K,

π_{K} (x_{i}) = \frac{1}{\sum_{h = 1}^{K} \exp (α_{h} + ϕ_{h} x_{i}^{T} β)} = \frac{1}{1 + \sum_{h = 1}^{K - 1} \exp (α_{h} + ϕ_{h} x_{i}^{T} β)} .

(10)

Note that these probabilities are equivalent to expressing the stereotype logit model using the Kth class as the baseline category; the j = 1, …, K − 1 logits are of the form

\log (\frac{π_{j} (x_{i})}{π_{K} (x_{i})}) = α_{j} + ϕ_{j} x_{i}^{T} β

(11)

such that the class-specific probabilities for classes j = 1, …, K − 1 are expressed as

π_{j} (x_{i}) = \exp (α_{j} + ϕ_{j} x_{j}^{T} β) π_{K} (x_{i}),

(12)

where

π_{K} (x_{i}) = \frac{1}{1 + \sum_{h = 1}^{K - 1} \exp (α_{h} + ϕ_{h} x_{i}^{T} β)} .

(13)

Substituting π_j(x_i) and π_K(x_i) into equation 2 and simplifying yields

\sum_{i = 1}^{n} (\sum_{k = 1}^{K - 1} y_{i k} (α_{k} - ϕ_{k} x_{i}^{T} β) + \log (1) - \log \sum_{h = 1}^{K} \exp (α_{h} - ϕ_{h} x_{i}^{T} β)) .

(14)

Because log(1) = 0, our final log-likelihood is

l (α, ϕ, β | x, y) = \sum_{i = 1}^{n} (\sum_{k = 1}^{K - 1} y_{i k} (α_{k} - ϕ_{k} x_{i}^{T} β) - \log \sum_{h = 1}^{K} \exp (α_{h} - ϕ_{h} x_{i}^{T} β))

(15)

or given the second form for π_K(x_i) in equation 10,

\log L (α, ϕ, β | x, y) = \sum_{i = 1}^{n} (\sum_{k = 1}^{K - 1} y_{i k} (α_{k} - ϕ_{k} x_{i}^{T} β) - \log (1 + \sum_{h = 1}^{K - 1} \exp (α_{h} - ϕ_{h} x_{i}^{T} β))) .

(16)

To ensure ordinality of the model,¹³ it was further required that ϕ₁ = 1, ϕ_K = 0, and ϕ₁ ≥ ϕ₂ ≥ ϕ₃≥ ··· ≥ ϕ_K-1 ≥ ϕ_K.

The first derivative of the log-likelihood with respect to β_p is given by

\begin{matrix} \frac{δ \log L}{δ β_{p}} = x_{p}^{T} (\sum_{h = 1}^{K - 1} y_{h} (ϕ_{h} - \frac{\sum_{h = 1}^{K - 1} ϕ_{h} \exp (α_{h} + ϕ_{h} x β)}{1 + \sum_{h = 1}^{K - 1} \exp (α_{h} + | ϕ_{h} x β)}) \\ - y_{K} \frac{\sum_{h = 1}^{K - 1} ϕ_{h} \exp (α_{h} + ϕ_{h} x β)}{1 + \sum_{h = 1}^{K - 1} \exp (α_{h} + ϕ_{h} x β)}) . \end{matrix}

(17)

The α terms are initialized as

α_{j} = \log (\frac{(\sum_{i = 1}^{n} y_{i j})}{\sum_{i = 1}^{n} y_{i K}}),

(18)

while ϕ are initialized as ϕ₁ = 1, ϕ₂, …, ϕ_K-1 = 0.1. If a no penalty subset is included, α, ϕ, and θ are estimated by maximum likelihood at the initial step and after each β update.

Implementation

The ordinalgmifs package was written in the R programing environment.¹⁵ The ordinal.gmifs function allows the user to specify a model formula, identify the matrix of covariates to be penalized in the model fitting algorithm using the x parameter, and additionally specify the model type (probability.model) and link function (link). The default is to fit a cumulative logit model, though allowable probability models include Cumulative, ForwardCR, BackwardCR, AdjCategory, and Stereotype while allowable links include logit, probit, and cloglog for the first three and loge and logit for the last two, respectively. The defaults for updating the penalized coefficients are epsilon = 0.001 and tol = 1e − 5. Our likelihood functions were written in R and tested by comparing our R output to output produced by the vglm R VGAM package for cumulative link, adjacent category, and forward and backward continuation ratio models and to STATA's slogit function and the rrvglm function in the R VGAM for the stereotype logit model using benchmark datasets for data where P < n.

Examples

The ordinalgmifs package includes example datasets having an ordinal response and a detailed tutorial as a vignette. The example data are a subset of subjects and CpG sites reported in the original paper where liver samples were assayed using the Illumina GoldenGate Methylation BeadArray Cancer Panel I.⁸ Technical replicate samples and matched cirrhotic samples from subjects with hepatocellular carcinoma (HCC) were removed to ensure all samples were independent. These data are in two formats: as a data.frame (hccframe) and as a BioConductor ExpressionSet (hccmethyl). However, for the demonstration provided here, we illustrate downloading the full dataset, GSE18081, from Gene Expression Omnibus, filtering on relevant criterion, and fitting a cumulative logit model to the independent subjects whose liver was either normal (n = 20) < cirrhotic but not having HCC (n = 16, cirrhosis non-HCC) < HCC (n = 20, HCC) using the CpG site methylation values as predictor variables. To download the GSE18081 dataset, we load the GEOquery R package before using the getGEO function.

> library(“GEOquery”)

> data<-getGEO(“GSE18081”) [[1]]

We then removed paired cirrhotic samples from patients that also contributed to tumor samples and replicate hybridizations so that only independent subjects are included in our dataset.

> data<-data [,ifelse (pData (data) $characteristics_ ch1==“disease state: Cirrhosis”, FALSE, TRUE)]

> data<-Data [,-grep (“Replicate”,pData (data) $title)]

> hccCancerPanel<-Data.frame(Tissue = factor(ifelse (pData (data) $characteristics_ch1==“disease state: Cirrhosis non-HCC”, “Cirrhosis non-HCC”,

ifelse (pData(data) $characteristics_ch1==“disease state: Normal”, “Normal”,

“HCC”)), ordered=TRUE, levels=c(“Normal”, “Cirrhosis non—HCC”, “HCC”)), t(exprs(data))) > rownames(hccCancerPanel)<-rownames (pData(data))

The class to be predicted (normal < cirrhosis non-HCC < HCC) is stored in our hccCancerPanel data.frame as Tissue along with the 1,505 CpG site methylation values. However, prior to model fitting, NA values should be imputed or removed from the data.frame. Our first filtering step applied removed 10 CpG sites that had at least one missing value. Additionally, we removed any CpG site that had a variance of 0 (n = 26), leaving 1,469 CpG sites for our predictive model.

> filter<-c(0,apply(hccCancerPanel[,-1],2,function(x) sum (is.na(x))))

> hccCancerPanel<-hccCancerPanel [,filter==0]

> filter2<-c(1,apply(hccCancerPanel[,-1],2,function(x) sd(x)))

> hccCancerPanel<-hccCancerPanel [, (1: dim (hccCancerPanel) [2]) [filter2!=0]]

To fit a model where all predictors are penalized, the model formula is specified to fit an intercept only model and the predictors to be penalized are specified using the x parameter. When fitting a penalized model, it is expected that more than one variable is included in the x parameter. The x parameter can either be a vector naming columns in the data.frame specified by the data parameter or be a data.frame name with the columns to include (or exclude) indicated by their (negative) index. By default, a cumulative logit model is fit when neither probability.model nor link is specified by the user. Because Tissue is the first variable in hccCancerPanel, we fit a model penalizing all CpG sites by specifying x = hccCancerPanel[,-1], which simply removes our ordinal outcome.

> library(“ordinalgmifs”)

> cumulative.logit<-ordinal.gmifs(Tissue ∼ 1, x = hccCancerPanel[,-1], data = hccCancerPanel)

Because the GMIFS procedure is incremental, the user may want to specify verbose = TRUE to print the step number in order to monitor the status of the model fitting procedure.

Methods including coef, plot, predict, fitted, print, and summary can be applied to ordinalgmifs model objects. Because the returned list differs depending on whether a no penalty subset is included or a stereotype logit model is fit, the print function returns the object names of the fitted object.

By default coef, predict, and summary extracts the relevant information from the step in the solution path that attained the minimum AIC. Because most of the features will have a coefficient estimate of zero, it is more convenient to extract the coefficient estimates and then examine those features with non-zero estimates. Here we see that in the AIC selected model there are 15 CpG sites with non-zero coefficient estimates.

Although the AIC is the default model selected, any step along the solution path can be extracted by specifying the step using the model.select parameter for these three functions. For example, the model attaining the minimum BIC can be extracted using

summary(cumulative.logit, model.select = which.min (cumulative.logit$BIC)).

Alternatively, the 150th step can be extracted using summary (cumulative.logit, model.select=150). Note that the α_j thresholds are labeled as (Intercept): 1, …, (Intercept):K − 1.

When examining the non-zero coefficient estimates, SOX17 and DDIT3 had the largest absolute values. Boxplots of the β values by tissue type revealed a monotonic relationship for both (Fig. 1).

Figure 1.

Boxplots of β values for SOX17 (left panel) and DDIT3 (right panel) by tissue type.

The plot function plots the solution path of the model fit. The vertical axis can be changed using the type parameter with allowable selections being trace (default), AIC, BIC, or logLik. Although there are default x-axis, y-axis, and titles provided for each plot, the user can modify these by supplying their own arguments to xlab, ylab, and main, respectively.

> plot(cumulative.logit)

The predict function (or equivalently, fitted) returns a list containing predicted, a matrix of the class probabilities from the fitted model, and class, the class having the maximum predicted probability from the fitted model. As with coef and summary, the predict function by default extracts the model that attained the minimum AIC, but predictions for any step along the solution path can be obtained by specifying the step using the model.select parameter.

For the AIC selected model, there were no misclassification errors. However, a fair way of estimating generalization error should be applied. When there are small sample sizes in one or more groups, K-fold cross-validation (CV) methods may not perform well as a means to estimate generalization error because of the random inclusion of samples into each of the folds. That is, multiple folds may include few if any subjects from the small classes. Therefore, here we have demonstrated N-fold CV for this dataset. Note that we include the drop = FALSE argument to preserve the dimension format of the object when only one subject comprises the test set. Note that we used the foreach and doSNOW packages for parallel processing to speed up our computations.

There were 8 of 56 misclassifications from the N-fold CV procedure, which yields a generalized misclassification rate of 14.3%.

Aside from a logit link, a probit or complementary log-log link can be used in conjunction with the cumulative link probability model. These three links are also available for probabilitymodel = “ForwardCR” and probability.model = “BackwardCR”. A stereotype logit model only uses a logit link, while an adjacent category model only uses a log_e link. Misspecifying the link for either a stereotype logit or adjacent category yields a warning that is printed to the R console, but only the correct link is used in the model fit.

Summary

Herein we presented an extension of the GMIFS method for predicting an ordinal response in high-dimensional covariates spaces. If interest lies in predicting an ordinal response and the number of covariates is small relative to the available sample size, we recommend using the VGAM package. However, our ordinalgmifs R package is capable of fitting models in high-dimensional covariate spaces and the penalization process better handles multicollinearity problems. We presented the flexibility of the GMIFS method by adapting it to fit cumulative link model, forward and backward continuation ratio models using either a logit, probit, or complementary log-log link, as well as the adjacent category and stereotype logit models. Functions for extracting coefficients, obtaining fitted probabilities, predicting class, plotting, and summarizing the fitted model object are also provided. Our ordinalgmifs package should be helpful when predicting an ordinal response for datasets where the number of covariates exceeds the number of available samples.

Author Contributions

Conceived and designed the methods: KJA, JH. Wrote the Software code: KJA, JH, QZ, KF, JGL, AEG. Wrote the first draft of the manuscript: KJA, JH. Contributed to the writing of the manuscript: QZ, KF, AEG. Agree with manuscript results and conclusions: KJA, JH, QZ, KF, AEG. Jointly developed the structure and arguments for the paper: KJA, JH. Made critical revisions and approved final version: KJA, JH, KF. All authors reviewed and approved of the final manuscript.

References

Park

M.Y.

, Hastie

L1-regularization path algorithm for generalized linear models.

J R Stat Soc B. 2007; 64: 659–77.

Friedman

J.H.

, Hastie

, Tibshirani

Regularization paths for generalized linear models via coordinate descent.

J Stat Softw. 2010; 33: 1–22.

Archer

K.J.

, Williams

A.A.A.

L1 penalized continuation ratio models for ordinal response prediction using high-dimensional datasets.

Stat Med. 2012; 31: 1464–74.

Tibshirani

Regression shrinkage and selection via the Lasso.

J R Stat Soc B. 1996; 58: 267–88.

Tibshirani

The lasso method for variable selection in the Cox model.

Stat Med. 1997; 16: 385–95.

Zou

, Hastie

Regularization and variable selection via the elastic net.

J R Stat Soc B. 2005; 67: 301–20.

Hastie

, Taylor

, Tibshirani

, Walther

Forward stagewise regression and the monotone lasso.

Electron J Stat. 2007; 1: 1–29.

Archer Kellie

, Mas Valeria

, Maluf Daniel

, Fisher Robert

High-throughput assessment of CpG site methylation for distinguishing between HCV-cirrhosis and HCV-associated hepatocellular carcinoma.

Mol Genet Genomics. 2010; 283: 341–9.

Efron

, Hastie

, Johnstone

, Tibshirani

Least angle regression.

Ann Stat. 2004; 32: 407–99.

10.

McCullagh

Regression models for ordinal data.

J R Stat Soc B. 1980; 42: 109–42.

11.

Agresti

Analysis of Ordinal Categorical Data. Hoboken, NJ: John Wiley & Sons; 2010.

12.

Bender

, Benner

Calculating ordinal regression models in SAS and S-Plus.

Biom J. 2000; 42: 677–99.

13.

Anderson

J.A.

Regression and ordered categorical variables.

JR Stat Soc Ser B. 1984; 46: 1–30.

14.

Greenland

Alternative models for ordinal logistic regression.

Stat Med. 1994; 13: 1665–77.

15.

R Core Team. R: A Language and Environment for Statistical Computing. Austria: R Foundation for Statistical Computing Vienna; 2013.