Nothing
#' Uses K-fold cross validation to obtain out-of-sample log-likelihood and
#' misclassification rates. Lambda is tuned within each cross validation fold.
#'
#' The data is divided into K folds. \code{ordinalNet} is fit \eqn{K} times, each time
#' leaving out one fold as a test set. For each of the \eqn{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.
#'
#' @param x Covariate matrix.
#' @param 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.
#' @param lambdaVals An optional user-specified lambda sequence (vector). If \code{NULL},
#' a sequence will be generated using the model fit to the full training data.
#' This default sequence is based on \code{nLambda} and \code{lambdaMinRatio},
#' which can be passed as additional arguments (otherwise \code{ordinalNet} default
#' values are used). 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 \code{lambdaMinRatio}.
#' @param folds An optional list, where each element is a vector of row indices
#' corresponding to a different cross validation fold. Indices correspond to rows
#' of the \code{x} matrix. Each index number should be used in exactly one fold.
#' If \code{NULL}, the data will be randomly divided into equally-sized partitions.
#' It is recommended to call \code{set.seed} before calling \code{ordinalNetCV}
#' for reproducibility.
#' @param nFolds Numer of cross validation folds. Only used if \code{folds=NULL}.
#' @param nFoldsCV Number of cross validation folds used to tune lambda for each
#' training set (i.e. within each training fold). Only used of \code{tuneMethod} is
#' "cvLoglik", "cvMisclass", "cvBrier", or "cvDevPct.
#' @param tuneMethod Method used to tune lambda for each training set (ie. within
#' each training fold). The "cvLoglik", "cvMisclass", "cvBrier", and "cvDevPct"
#' methods use cross validation with \code{nFoldsCV} folds and select the
#' lambda value with the best average out-of-sample performance. The "aic" and "bic"
#' methods are less computationally intensive because they do not require the
#' model to be fit multiple times.
#' Note that for the methods that require cross validation, the fold splits are
#' determined randomly and cannot be specified by the user. The \code{set.seed()}
#' function should be called prior to \code{ordinalNetCV} for reproducibility.
#' @param printProgress Logical. If \code{TRUE} the fitting progress is printed
#' to the terminal.
#' @param warn Logical. If \code{TRUE}, the following warning message is displayed
#' when fitting a cumulative probability model with \code{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.
#' @param ... Other arguments (besides \code{x}, \code{y}, \code{lambdaVals}, and \code{warn})
#' passed to \code{ordinalNet}.
#'
#' @details
#' \itemize{
#' \item The fold partition splits can be passed by the user via the \code{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 \code{set.seed}
#' function should be called prior to \code{ordinalNetCV} for reproducibility.
#' \item A sequence of lambda values can be passed by the user via the
#' \code{lambdaVals} argument. By default, the sequence is generated by first
#' fitting the model to the full data set (this sequence is determined by the
#' \code{nLambda} and \code{lambdaMinRatio} arguments of \code{ordinalNet}).
#' \item The \code{standardize} argument of \code{ordinalNet} can be modified through
#' the additional arguments (...). If \code{standardize=TRUE}, then the data are scaled
#' within each cross validation fold. If \code{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.
#' }
#'
#' @return
#' An S3 object of class "ordinalNetCV", which contains the following:
#' \describe{
#' \item{loglik}{Vector of out-of-sample log-likelihood values.
#' Each value corresponds to a different fold.}
#' \item{misclass}{Vector of out-of-sample misclassificaton rates.
#' Each value corresponds to a different fold.}
#' \item{brier}{Vector of out-of-sample Brier scores.
#' Each value corresponds to a different fold.}
#' \item{devPct}{Vector of out-of-sample percentages of deviance explained.
#' Each value corresponds to a different fold.}
#' \item{bestLambdaIndex}{The index of the value within the lambda sequence
#' selected for each fold by the tuning method.}
#' \item{lambdaVals}{The sequence of lambda values used for all cross validation folds.}
#' \item{folds}{A list containing the index numbers of each fold.}
#' \item{fit}{An object of class "ordinalNet", resulting from fitting
#' \code{ordinalNet} to the entire dataset.}
#' }
#'
#' @examples
#' \dontrun{
#' # 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)
#' }
#'
#' @export
ordinalNetCV <- function(x, y, lambdaVals=NULL, folds=NULL, nFolds=5, nFoldsCV=5,
tuneMethod=c("cvLoglik", "cvMisclass", "cvBrier", "cvDevPct",
"aic", "bic"),
printProgress=TRUE, warn=TRUE, ...)
{
tuneMethod <- match.arg(tuneMethod)
# cvID := indicator to use cross validation within folds
cvID <- tuneMethod %in% c("cvLoglik", "cvMisclass", "cvBrier", "cvDevPct")
# tuneMethod := argument passed to ordinalNetTune
if (tuneMethod == "cvLoglik") cvCriterion <- "loglik"
if (tuneMethod == "cvMisclass") cvCriterion <- "misclass"
if (tuneMethod == "cvBrier") cvCriterion <- "brier"
if (tuneMethod == "cvDevPct") cvCriterion <- "devPct"
# Argument checks
if (is.matrix(y) && any(rowSums(y)!=1))
warning(paste0("Data is split by row for cross validation, but note that ",
"y matrix rows have different weights. Be sure this is what you want."))
if (!is.null(folds) && length(folds)<2)
stop(paste0("\'folds\' should be a list of at least two vectors. ",
"Each vector should contain indices of a cross validation fold. ",
"Each index from 1:nrow(x) should be used exactly once."))
if (!is.null(folds) && !setequal(unlist(folds), 1:nrow(x)))
stop("\'folds\' should include each index from 1:nrow(x) exactly once.")
yMat <- if (is.matrix(y)) y else yFactorToMatrix(y) # for computing log-likelihood
if (printProgress) cat("Fitting ordinalNet on full training data\n")
fit <- ordinalNet(x, y, lambdaVals=lambdaVals, warn=warn, ...)
if (is.null(lambdaVals)) lambdaVals <- fit$lambdaVals
if (is.null(folds))
{
n <- nrow(x)
randIndex <- sample(n)
folds <- split(randIndex, rep(1:nFolds, length.out=n))
} else
{
nFolds <- length(folds)
}
nLambda <- length(lambdaVals)
loglik <- misclass <- brier <- devPct <- bestLambdaIndex <- rep(NA, nFolds)
names(loglik) <- names(misclass) <- names(brier) <- names(devPct) <-
names(bestLambdaIndex) <- paste0("fold", 1:nFolds)
for (i in 1:nFolds)
{
testFold <- folds[[i]]
xTrain <- x[-testFold, , drop=FALSE]
xTest <- x[testFold, , drop=FALSE]
yTrain <- if (is.matrix(y)) y[-testFold, , drop=FALSE] else y[-testFold]
yMatTest <- yMat[testFold, , drop=FALSE]
if (printProgress) cat("Fitting ordinalNet on fold", i, "of", nFolds, '\n')
if (cvID)
{
fitTrainCV <- ordinalNetTune(xTrain, yTrain, lambdaVals=lambdaVals, folds=NULL,
nFolds=5, printProgress=FALSE, warn=FALSE, ...)
fitTrain <- fitTrainCV$fit
if (cvCriterion %in% c("loglik", "devPct"))
wm <- which.max
if (cvCriterion %in% c("misclass", "brier"))
wm <- which.min
bestLambdaIndex[[i]] <- wm(rowMeans(fitTrainCV[[cvCriterion]]))
} else # tuneMethod is either "aic" or "bic"
{
fitTrain <- ordinalNet(xTrain, yTrain, lambdaVals=lambdaVals, warn=FALSE, ...)
bestLambdaIndex[[i]] <- which.min(fitTrain[[tuneMethod]])
}
pHatFull <- predict.ordinalNet(fitTrain, newx=xTest, type="response", whichLambda=bestLambdaIndex[[i]])
pHat <- pHatFull[, -ncol(pHatFull), drop=FALSE]
loglik[i] <- getLoglik(pHat, yMatTest)
misclass[i] <- getMisclass(pHat, yMatTest)
brier[i] <- getBrier(pHat, yMatTest)
loglikNull <- getLoglikNull(yMatTest)
devPct[i] <- 1 - loglik[i] / loglikNull
}
if (printProgress) cat("Done\n")
out <- list(loglik=loglik, misclass=misclass, brier=brier, devPct=devPct,
bestLambdaIndex=bestLambdaIndex, lambdaVals=lambdaVals, folds=folds, fit=fit)
class(out) <- "ordinalNetCV"
out
}
#' Summary method for an "ordinalNetCV" object.
#'
#' 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.
#'
#' @param object An "ordinalNetCV" S3 object
#' @param ... Not used. Additional summary arguments.
#'
#' @return 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.
#'
#' @seealso
#' \code{\link{ordinalNetCV}}
#'
#' @examples
#' # See ordinalNetCV() documentation for examples.
#'
#' @export
summary.ordinalNetCV <- function(object, ...)
{
lambda <- object$lambdaVals[object$bestLambdaIndex]
loglik <- object$loglik
misclass <- object$misclass
brier <- object$brier
devPct <- object$devPct
data.frame(lambda=lambda, loglik=loglik, misclass=misclass, brier=brier, devPct=devPct)
}
#' Print method for an "ordinalNetCV" object.
#'
#' Prints the data frame returned by the \code{summary.ordinalNetCV()} method.
#'
#' @param x An "ordinalNetCV" S3 object
#' @param ... Not used. Additional print arguments.
#'
#' @seealso
#' \code{\link{ordinalNetCV}}
#'
#' @examples
#' # See ordinalNetCV() documentation for examples.
#'
#' @export
print.ordinalNetCV <- function(x, ...)
{
cat("\nCross validation summary:\n\n")
print(summary.ordinalNetCV(x))
cat("\n")
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.