Method to extract fitted coefficients from an "ordinalNet" object.

Description

Method to extract fitted coefficients from an "ordinalNet" object.

Usage

## S3 method for class 'ordinalNet'
coef(
  object,
  matrix = FALSE,
  whichLambda = NULL,
  criteria = c("aic", "bic"),
  ...
)

Arguments

Value

The object returned depends on matrix.

See Also

ordinalNet

Examples

# See ordinalNet() documentation for examples.

Ordinal regression models with elastic net penalty

Description

Fits ordinal regression models with elastic net penalty by coordinate descent. Supported model families include cumulative probability, stopping ratio, continuation ratio, and adjacent category. These families are a subset of vector glm's which belong to a model class we call the elementwise link multinomial-ordinal (ELMO) class. Each family in this class links a vector of covariates to a vector of class probabilities. Each of these families has a parallel form, which is appropriate for ordinal response data, as well as a nonparallel form that is appropriate for an unordered categorical response, or as a more flexible model for ordinal data. The parallel model has a single set of coefficients, whereas the nonparallel model has a set of coefficients for each response category except the baseline category. It is also possible to fit a model with both parallel and nonparallel terms, which we call the semi-parallel model. The semi-parallel model has the flexibility of the nonparallel model, but the elastic net penalty shrinks it toward the parallel model.

Usage

ordinalNet(
  x,
  y,
  alpha = 1,
  standardize = TRUE,
  penaltyFactors = NULL,
  positiveID = NULL,
  family = c("cumulative", "sratio", "cratio", "acat"),
  reverse = FALSE,
  link = c("logit", "probit", "cloglog", "cauchit"),
  customLink = NULL,
  parallelTerms = TRUE,
  nonparallelTerms = FALSE,
  parallelPenaltyFactor = 1,
  lambdaVals = NULL,
  nLambda = 20,
  lambdaMinRatio = 0.01,
  includeLambda0 = FALSE,
  alphaMin = 0.01,
  pMin = 1e-08,
  stopThresh = 1e-08,
  threshOut = 1e-08,
  threshIn = 1e-08,
  maxiterOut = 100,
  maxiterIn = 100,
  printIter = FALSE,
  printBeta = FALSE,
  warn = TRUE,
  keepTrainingData = TRUE
)

Arguments

Details

The ordinalNet function fits regression models for a categorical response variable with K+1 levels. Conditional on the covariate vector x_i (the i^{th} row of x), each observation has a vector of K+1 class probabilities (p_{i1}, \ldots, p_{i(K+1)}). These probabilities sum to one, and can therefore be parametrized by p_i = (p_{i1}, \ldots, p_{iK}). The probabilities are mapped to a set of K quantities \delta_i = (\delta_{i1}, \ldots, \delta_{iK}) \in (0, 1)^K, which depends on the choice of model family. The elementwise link function maps \delta_i to a set of K linear predictors. Together, the family and link specifiy a link function between p_i and \eta_i.

Model families:

Let Y denote the random response variable for a single observation, conditional on the covariates values of the observation. The random variable Y is discrete with support {1, \ldots, K+1}. The following model families are defined according to these mappings between the class probabilities and the values \delta_1, \ldots, \delta_K:

Cumulative probability

\delta_j = P(Y \le j)

Reverse cumulative probability

\delta_j = P(Y \ge j + 1)

Stopping ratio

\delta_j = P(Y = j | Y \ge j)

Reverse stopping ratio

\delta_j = P(Y=j + 1 | Y \le j + 1)

Continuation ratio

\delta_j = P(Y > j | Y \ge j)

Reverse continuation ratio

\delta_j = P(Y < j | Y \le j)

Adjacent category

\delta_j = P(Y = j + 1 | j \le Y \le j+1)

Reverse adjacent category

\delta_j = P(Y = j | j \le Y \le j+1)

Parallel, nonparallel, and semi-parallel model forms:

Models within each of these families can take one of three forms, which have different definitions for the linear predictor \eta_i. Suppose each x_i has length P. Let b be a length P vector of regression coefficients. Let B be a P \times K matrix of regression coefficient. Let b_0 be a vector of K intercept terms. The three model forms are the following:

Parallel

\eta_i = b_0 + b^T x_i (parallelTerms=TRUE, nonparallelTerms=FALSE)

Nonparallel

\eta_i = b_0 + B^T x_i (parallelTerms=FALSE, nonparallelTerms=TRUE)

Semi-parallel

\eta_i = b_0 + b^T x_i + B^T x_i (parallelTerms=TRUE, nonparallelTerms=TRUE)

The parallel form has the defining property of ordinal models, which is that a single linear combination b^T x_i shifts the cumulative class probabilities P(Y \le j) in favor of either higher or lower categories. The linear predictors are parallel because they only differ by their intercepts (b_0). The nonparallel form is a more flexible model, and it does not shift the cumulative probabilities together. The semi-parallel model is equivalent to the nonparallel model, but the elastic net penalty shrinks the semi-parallel coefficients toward a common value (i.e. the parallel model), as well as shrinking all coefficients toward zero. The nonparallel model, on the other hand, simply shrinks all coefficients toward zero. When the response categories are ordinal, any of the three model forms could be applied. When the response categories are unordered, only the nonparallel model is appropriate.

Elastic net penalty:

The elastic net penalty is defined for each model form as follows. \lambda and \alpha are the usual elastic net tuning parameters, where \lambda determines the degree to which coefficients are shrunk toward zero, and \alpha specifies the amound of weight given to the L1 norm and squared L2 norm penalties. Each covariate is allowed a unique penalty factor c_j, which is specified with the penaltyFactors argument. By default c_j = 1 for all j. The semi-parallel model has a tuning parameter \rho which determines the degree to which the parallel coefficients are penalized. Small values of \rho will result in a fit closer to the parallel model, and large values of \rho will result in a fit closer to the nonparallel model.

Parallel

\lambda \sum_{j=1}^P c_j \{ \alpha |b_j| + \frac{1}{2} (1-\alpha) b_j^2 \}

Nonparallel

\lambda \sum_{j=1}^P c_j \{ \sum_{k=1}^K \alpha |B_{jk}| + \frac{1}{2} (1-\alpha) B_{jk}^2 \}

Semi-parallel

\lambda [ \rho \sum_{j=1}^P c_j \{ \alpha |b_j| + \frac{1}{2} (1-\alpha) b_j^2 \} + \sum_{j=1}^P c_j \{ \sum_{k=1}^K \alpha |B_{jk}| + \frac{1}{2} (1-\alpha) B_{jk}^2 \}]

ordinalNet minimizes the following objective function. Let N be the number of observations, which is defined as the sum of the y elements when y is a matrix.

objective = -1/N*loglik + penalty

Value

An object with S3 class "ordinalNet". Model fit information can be accessed through the coef, predict, and summary methods.

coefs

Matrix of coefficient estimates, with each row corresponding to a lambda value. (If covariates were scaled with standardize=TRUE, the coefficients are returned on the original scale).

lambdaVals

Sequence of lambda values. If user passed a sequence to the lambdaVals, then it is this sequence. If lambdaVals argument was NULL, then it is the sequence generated.

loglik

Log-likelihood of each model fit.

nNonzero

Number of nonzero coefficients of each model fit, including intercepts.

aic

AIC, defined as -2*loglik + 2*nNonzero.

bic

BIC, defined as -2*loglik + log(N)*nNonzero.

devPct

Percentage deviance explained, defined as 1 - loglik/loglik_0, where loglik_0 is the log-likelihood of the null model.

iterOut

Number of coordinate descent outer loop iterations until convergence for each lambda value.

iterIn

Number of coordinate descent inner loop iterations on last outer loop for each lambda value.

dif

Relative improvement in objective function on last outer loop for each lambda value. Can be used to diagnose convergence issues. If iterOut reached maxiterOut and dif is large, then maxiterOut should be increased. If dif is negative, this means the objective did not improve between successive iterations. This usually only occurs when the model is saturated and/or close to convergence, so a small negative value is not of concern. (When this happens, the algorithm is terminated for the current lambda value, and the coefficient estimates from the previous outer loop iteration are returned.)

nLev

Number of response categories.

nVar

Number of covariates in x.

xNames

Covariate names.

args

List of arguments passed to the ordinalNet function.

Examples

# Simulate x as independent standard normal
# Simulate y|x from a parallel cumulative logit (proportional odds) model
set.seed(1)
n <- 50
intercepts <- c(-1, 1)
beta <- c(1, 1, 0, 0, 0)
ncat <- length(intercepts) + 1  # number of response categories
p <- length(beta)  # number of covariates
x <- matrix(rnorm(n*p), ncol=p)  # n x p covariate matrix
eta <- c(x %*% beta) + matrix(intercepts, nrow=n, ncol=ncat-1, byrow=TRUE)
invlogit <- function(x) 1 / (1+exp(-x))
cumprob <- t(apply(eta, 1, invlogit))
prob <- cbind(cumprob, 1) - cbind(0, cumprob)
yint <- apply(prob, 1, function(p) sample(1:ncat, size=1, prob=p))
y <- as.factor(yint)

# Fit parallel cumulative logit model
fit1 <- ordinalNet(x, y, family="cumulative", link="logit",
                   parallelTerms=TRUE, nonparallelTerms=FALSE)
summary(fit1)
coef(fit1)
coef(fit1, matrix=TRUE)
predict(fit1, type="response")
predict(fit1, type="class")

# Fit nonparallel cumulative logit model
fit2 <- ordinalNet(x, y, family="cumulative", link="logit",
                   parallelTerms=FALSE, nonparallelTerms=TRUE)
fit2
coef(fit2)
coef(fit2, matrix=TRUE)
predict(fit2, type="response")
predict(fit2, type="class")

# Fit semi-parallel cumulative logit model (with both parallel and nonparallel terms)
fit3 <- ordinalNet(x, y, family="cumulative", link="logit",
                   parallelTerms=TRUE, nonparallelTerms=TRUE)
fit3
coef(fit3)
coef(fit3, matrix=TRUE)
predict(fit3, type="response")
predict(fit3, type="class")

Uses K-fold cross validation to obtain out-of-sample log-likelihood and misclassification rates. Lambda is tuned within each cross validation fold.

Description

The data is divided into K folds. ordinalNet is fit K times, each time leaving out one fold as a test set. For each of the K model fits, lambda can be tuned by AIC or BIC, or cross validation. If cross validation is used, the user can choose whether to user the best average out-of-sample log-likelihood, misclassification rate, Brier score, or percentage of deviance explained. The user can also choose the number of cross validation folds to use for tuning. Once the model is tuned, the out of sample log-likelihood, misclassification rate, Brier score, and percentage of deviance explained are calculated on the held out test set.

Usage

ordinalNetCV(
  x,
  y,
  lambdaVals = NULL,
  folds = NULL,
  nFolds = 5,
  nFoldsCV = 5,
  tuneMethod = c("cvLoglik", "cvMisclass", "cvBrier", "cvDevPct", "aic", "bic"),
  printProgress = TRUE,
  warn = TRUE,
  ...
)

Arguments

Details

• The fold partition splits can be passed by the user via the folds argument. By default, the data are randomly divided into equally-sized partitions. Note that if lambda is tuned by cross validation, the fold splits are determined randomly and cannot be specified by the user. The set.seed function should be called prior to ordinalNetCV for reproducibility.

• A sequence of lambda values can be passed by the user via the lambdaVals argument. By default, the sequence is generated by first fitting the model to the full data set (this sequence is determined by the nLambda and lambdaMinRatio arguments of ordinalNet).

• The standardize argument of ordinalNet can be modified through the additional arguments (...). If standardize=TRUE, then the data are scaled within each cross validation fold. If standardize=TRUE and lambda is tuned by cross validation, then the data are also scaled within each tuning sub-fold. This is done because scaling is part of the statistical procedure and should be repeated each time the procedure is applied.

Value

An S3 object of class "ordinalNetCV", which contains the following:

loglik

Vector of out-of-sample log-likelihood values. Each value corresponds to a different fold.

misclass

Vector of out-of-sample misclassificaton rates. Each value corresponds to a different fold.

brier

Vector of out-of-sample Brier scores. Each value corresponds to a different fold.

devPct

Vector of out-of-sample percentages of deviance explained. Each value corresponds to a different fold.

bestLambdaIndex

The index of the value within the lambda sequence selected for each fold by the tuning method.

lambdaVals

The sequence of lambda values used for all cross validation folds.

folds

A list containing the index numbers of each fold.

fit

An object of class "ordinalNet", resulting from fitting ordinalNet to the entire dataset.

Examples

## Not run: 
# Simulate x as independent standard normal
# Simulate y|x from a parallel cumulative logit (proportional odds) model
set.seed(1)
n <- 50
intercepts <- c(-1, 1)
beta <- c(1, 1, 0, 0, 0)
ncat <- length(intercepts) + 1  # number of response categories
p <- length(beta)  # number of covariates
x <- matrix(rnorm(n*p), ncol=p)  # n x p covariate matrix
eta <- c(x %*% beta) + matrix(intercepts, nrow=n, ncol=ncat-1, byrow=TRUE)
invlogit <- function(x) 1 / (1+exp(-x))
cumprob <- t(apply(eta, 1, invlogit))
prob <- cbind(cumprob, 1) - cbind(0, cumprob)
yint <- apply(prob, 1, function(p) sample(1:ncat, size=1, prob=p))
y <- as.factor(yint)

# Evaluate out-of-sample performance of the  cumulative logit model
# when lambda is tuned by cross validation (best average out-of-sample log-likelihood)
cv <- ordinalNetCV(x, y, tuneMethod="cvLoglik")
summary(cv)

## End(Not run)

Uses K-fold cross validation to obtain out-of-sample log-likelihood and misclassification rates for a sequence of lambda values.

Description

The data is divided into K folds. ordinalNet is fit K times (K=nFolds), each time leaving out one fold as a test set. The same sequence of lambda values is used each time. The out-of-sample log-likelihood, misclassification rate, Brier score, and percentage of deviance explained are obtained for each lambda value from the held out test set. It is up to the user to determine how to tune the model using this information.

Usage

ordinalNetTune(
  x,
  y,
  lambdaVals = NULL,
  folds = NULL,
  nFolds = 5,
  printProgress = TRUE,
  warn = TRUE,
  ...
)

Arguments

Details

• The fold partition splits can be passed by the user via the folds argument. By default, the data are randomly divided into equally-sized partitions. The set.seed function should be called prior to ordinalNetCV for reproducibility.

• A sequence of lambda values can be passed by the user via the lambdaVals argument. By default, the sequence is generated by first fitting the model to the full data set (this sequence is determined by the nLambda and lambdaMinRatio arguments of ordinalNet).

• The standardize argument of ordinalNet can be modified through the additional arguments (...). If standardize=TRUE, then the data are scaled within each cross validation fold. This is done because scaling is part of the statistical procedure and should be repeated each time the procedure is applied.

Value

An S3 object of class "ordinalNetTune", which contains the following:

loglik

Matrix of out-of-sample log-likelihood values. Each row corresponds to a lambda value, and each column corresponds to a fold.

misclass

Matrix of out-of-sample misclassificaton rates. Each row corresponds to a lambda value, and each column corresponds to a fold.

brier

Matrix of out-of-sample Brier scores. Each row corresponds to a lambda value, and each column corresponds to a fold.

devPct

Matrix of out-of-sample percentages of deviance explained. Each row corresponds to a lambda value, and each column corresponds to a fold.

lambdaVals

The sequence of lambda values used for all cross validation folds.

folds

A list containing the index numbers of each fold.

fit

An object of class "ordinalNet", resulting from fitting ordinalNet to the entire dataset.

Examples

## Not run: 
# Simulate x as independent standard normal
# Simulate y|x from a parallel cumulative logit (proportional odds) model
set.seed(1)
n <- 50
intercepts <- c(-1, 1)
beta <- c(1, 1, 0, 0, 0)
ncat <- length(intercepts) + 1  # number of response categories
p <- length(beta)  # number of covariates
x <- matrix(rnorm(n*p), ncol=p)  # n x p covariate matrix
eta <- c(x %*% beta) + matrix(intercepts, nrow=n, ncol=ncat-1, byrow=TRUE)
invlogit <- function(x) 1 / (1+exp(-x))
cumprob <- t(apply(eta, 1, invlogit))
prob <- cbind(cumprob, 1) - cbind(0, cumprob)
yint <- apply(prob, 1, function(p) sample(1:ncat, size=1, prob=p))
y <- as.factor(yint)

# Fit parallel cumulative logit model; select lambda by cross validation
tunefit <- ordinalNetTune(x, y)
summary(tunefit)
plot(tunefit)
bestLambdaIndex <- which.max(rowMeans(tunefit$loglik))
coef(tunefit$fit, whichLambda=bestLambdaIndex, matrix=TRUE)
predict(tunefit$fit, whichLambda=bestLambdaIndex)

## End(Not run)

Plot method for "ordinalNetTune" object.

Description

Plots the average out-of-sample log-likelihood, misclassification rate, Brier score, or percentage of deviance explained for each lambda value in the solution path. The averae is taken over all cross validation folds.

Usage

## S3 method for class 'ordinalNetTune'
plot(x, type = c("loglik", "misclass", "brier", "devPct"), ...)

Arguments

See Also

ordinalNetTune

Examples

# See ordinalNetTune() documentation for examples.

Predict method for an "ordinalNet" object

Description

Obtains predicted probabilities, predicted class, or linear predictors.

Usage

## S3 method for class 'ordinalNet'
predict(
  object,
  newx = NULL,
  whichLambda = NULL,
  criteria = c("aic", "bic"),
  type = c("response", "class", "link"),
  ...
)

Arguments

Value

The object returned depends on type.

See Also

ordinalNet

Examples

# See ordinalNet() documentation for examples.

Print method for an "ordinalNet" object.

Description

Prints the data frame returned by the summary.ordinalNet() method.

Usage

## S3 method for class 'ordinalNet'
print(x, ...)

Arguments

See Also

ordinalNet

Examples

# See ordinalNet() documentation for examples.

Print method for an "ordinalNetCV" object.

Description

Prints the data frame returned by the summary.ordinalNetCV() method.

Usage

## S3 method for class 'ordinalNetCV'
print(x, ...)

Arguments

See Also

ordinalNetCV

Examples

# See ordinalNetCV() documentation for examples.

Print method for an "ordinalNetTune" object.

Description

Prints the data frame returned by the summary.ordinalNetTune() method.

Usage

## S3 method for class 'ordinalNetTune'
print(x, ...)

Arguments

See Also

ordinalNetTune

Examples

# See ordinalNetTune() documentation for examples.

Summary method for an "ordinalNet" object.

Description

Provides a data frame which summarizes the model fit at each lambda value in the solution path.model fit summary as a data frame.

Usage

## S3 method for class 'ordinalNet'
summary(object, ...)

Arguments

Value

A data frame containing a record for each lambda value in the solution path. Each record contains the following fields: lambda value, degrees of freedom (number of nonzero parameters), log-likelihood, AIC, BIC, and percent deviance explained.

See Also

ordinalNet

Examples

# See ordinalNet() documentation for examples.

Summary method for an "ordinalNetCV" object.

Description

Provides a data frame which summarizes the cross validation results, which can be used as an estimate of the out-of-sample performance of a model tuned by a particular method.

Usage

## S3 method for class 'ordinalNetCV'
summary(object, ...)

Arguments

Value

A data frame containing a record for each cross validation fold. Each record contains the following: lambda value, log-likelihood, misclassification rate, Brier score, and percentage of deviance explained.

See Also

ordinalNetCV

Examples

# See ordinalNetCV() documentation for examples.

Summary method for an "ordinalNetTune" object.

Description

Provides a data frame which summarizes the cross validation results and may be useful for selecting an appropriate value for the tuning parameter lambda.

Usage

## S3 method for class 'ordinalNetTune'
summary(object, ...)

Arguments

Value

A data frame containing a record for each lambda value in the solution path. Each record contains the following: lambda value, average log-likelihood, average misclassification rate, average Brier score, and average percentage of deviance explained. Averages are taken across all cross validation folds.

See Also

ordinalNetTune

Examples

# See ordinalNetTune() documentation for examples.