diff --git a/R/classes.R b/R/classes.R index 4f96869..64f2906 100644 --- a/R/classes.R +++ b/R/classes.R @@ -13,9 +13,6 @@ setOldClass("PrecisionPathways") # Union of A Function and NULL setClassUnion("functionOrNULL", c("function", "NULL")) -# Union of a Function and a List of Functions. Useful for allowing ensemble feature selection. -setClassUnion("functionOrChraracterOrList", c("function", "character", "list")) - # Union of A Numeric Value and NULL setClassUnion("numericOrNULL", c("numeric", "NULL")) @@ -240,20 +237,30 @@ setMethod("show", "TransformParams", function(object) { cat("An object of class 'TransformParams'.\n") - index <- na.omit(match("Transform Name", object@characteristics[, "characteristic"])) - transText <- object@characteristics[index, "characteristic"] - cat("Transform Name: ", object@characteristics[index, "value"], ".\n", sep = '') - - otherInfo <- object@characteristics[-index, ] - if(nrow(otherInfo) > 0) - { - for(rowIndex in 1:nrow(otherInfo)) - { - cat(otherInfo[rowIndex, "characteristic"], ": ", otherInfo[rowIndex, "value"], ".\n", sep = '') - } - } + .showCharacteristics(object@characteristics, "Transform Name") }) +# Prints the characteristics named by IDs first, then any others, one per line. +.showCharacteristics <- function(characteristics, IDs) +{ + if(ncol(characteristics) == 0 || nrow(characteristics) == 0) return(invisible(NULL)) + index <- match(IDs, characteristics[, "characteristic"]) + index <- c(index[!is.na(index)], setdiff(seq_len(nrow(characteristics)), index)) + for(rowIndex in index) + cat(characteristics[rowIndex, "characteristic"], ": ", characteristics[rowIndex, "value"], ".\n", sep = '') + invisible(NULL) +} + +# The name of a feature selection, training or prediction function, as shown in characteristics. +# Functions of the package have a registered name. Other functions are described generically. +.functionDisplayName <- function(aFunction, default) +{ + functionName <- attr(aFunction, "name") + if(is.null(functionName)) return(default) + displayName <- .ClassifyRenvir[["functionsTable"]][.ClassifyRenvir[["functionsTable"]][, "character"] == functionName, "name"] + if(length(displayName) == 0) functionName else displayName[1] +} + ##### FeatureSetCollection ##### #' @docType class #' @exportClass FeatureSetCollection @@ -425,7 +432,7 @@ setMethod("show", "FeatureSetCollection", }) setsText <- paste(names(object@sets)[minIndex:lastIndex], featuresConcatenated, sep = ": ") setsText <- paste(setsText, collapse = '\n') - cat(setsText) + cat(setsText, "\n", sep = '') } } ) @@ -450,6 +457,9 @@ setClassUnion("FeatureSetCollectionOrNULL", c("FeatureSetCollection", "NULL")) setClassUnion("functionOrList", c("function", "list")) setClassUnion("characterOrList", c("character", "list")) +# Union of a function, a character keyword and a list of either. Allows ensemble feature selection. +setClassUnion("functionOrCharacterOrList", c("function", "character", "list")) + ##### SelectParams ##### #' @exportClass SelectParams @@ -477,13 +487,13 @@ setClassUnion("SelectParamsOrNULL", c("SelectParams", "NULL")) #' @name SelectParams #' @rdname SelectParams-class #' @aliases SelectParams SelectParams-class SelectParams,missing-method -#' SelectParams,characterOrList-method +#' SelectParams,functionOrCharacterOrList-method #' @docType class #' @section Constructor: #' \describe{ #' \item{\code{SelectParams(featureRanking, characteristics = DataFrame(), nFeatures = 20, minPresence = 1, intermediate = character(0), subsetToSelections = TRUE, tuneParams = list(nFeatures = seq(10, 100, 10)), ...)}}{Creates a \code{SelectParams} object which stores the function(s) which will do the selection and parameters that the function will use.\cr #' \describe{ -#' \item{\code{featureRanking}}{A character keyword referring to a registered feature ranking function. See \code{\link{available}} for valid keywords.} +#' \item{\code{featureRanking}}{A character keyword referring to a registered feature ranking function, a feature ranking function or a list of either. See \code{\link{available}} for valid keywords. A ranking function must accept the training measurements and outcome as its first two arguments and return the features ranked from best to worst. The keyword \code{"none"} returns \code{NULL}, which \code{\link{ModellingParams}} accepts as no feature selection.} #' \item{\code{characteristics}}{A \code{\link{DataFrame}} describing the characteristics of feature selection to be done. First column must be named \code{"charateristic"} and second column must be named \code{"value"}. If using wrapper functions for feature selection in this package, the feature selection name will automatically be generated and therefore it is not necessary to specify it.} #' \item{\code{nFeatures}}{Default: \code{20}. The number of top-ranked features to choose. Can also be \code{NULL} if a vector of top numbers is specified to \code{tuneParams} for the list element named \code{nFeatures}.} #' \item{\code{minPresence}}{Default: \code{1}. If a list of functions was provided, how many of those must a feature have been selected by to be used in classification. 1 is equivalent to a set union and a number the same length as \code{featureSelection} is equivalent to set intersection.} @@ -517,19 +527,32 @@ standardGeneric("SelectParams")) #' @rdname SelectParams-class #' @usage NULL #' @export -setMethod("SelectParams", c("characterOrList"), +setMethod("SelectParams", c("functionOrCharacterOrList"), function(featureRanking, characteristics = DataFrame(), nFeatures = 20, minPresence = 1, intermediate = character(0), subsetToSelections = TRUE, tuneParams = NULL, ...) { - if(is.character(featureRanking)) featureRanking <- .selectionKeywordToFunction(featureRanking) else featureRanking <- lapply(featureRanking, .selectionKeywordToFunction) + if(is.character(featureRanking) && length(featureRanking) == 1 && featureRanking == "none") + return(NULL) # No feature selection. + toFunction <- function(ranking) + { + if(is.function(ranking)) + { # Feature selection identifies functions by their name. + if(is.null(attr(ranking, "name"))) attr(ranking, "name") <- "User-specified Ranking" + return(ranking) + } + rankingFunction <- .selectionKeywordToFunction(ranking) + if(is.null(rankingFunction)) + stop("'", ranking, "' is not a feature selection keyword. See available(\"selectionMethod\") for valid keywords.") + rankingFunction + } + if(is.list(featureRanking)) featureRanking <- lapply(featureRanking, toFunction) else featureRanking <- toFunction(featureRanking) if(!is.list(featureRanking) && (ncol(characteristics) == 0 || !"Selection Name" %in% characteristics[, "characteristic"])) { - if(!is.null(featureRanking)) - characteristics <- rbind(characteristics, S4Vectors::DataFrame(characteristic = "Selection Name", value = .ClassifyRenvir[["functionsTable"]][.ClassifyRenvir[["functionsTable"]][, "character"] == attr(featureRanking, "name"), "name"])) + characteristics <- rbind(characteristics, S4Vectors::DataFrame(characteristic = "Selection Name", value = .functionDisplayName(featureRanking, "User-specified Ranking"))) } if(is.list(featureRanking) && (ncol(characteristics) == 0 || !"Ensemble Selection" %in% characteristics[, "characteristic"])) { - selectMethodNames <- unlist(lapply(featureRanking, function(rankingFunction) .ClassifyRenvir[["functionsTable"]][.ClassifyRenvir[["functionsTable"]][, "character"] == attr(rankingFunction, "name"), "name"])) + selectMethodNames <- unlist(lapply(featureRanking, .functionDisplayName, "User-specified Ranking")) characteristics <- rbind(characteristics, S4Vectors::DataFrame(characteristic = "Ensemble Selection", value = paste(selectMethodNames, collapse = ", "))) } if(!is.null(tuneParams[["nFeatures"]])) nFeatures <- NULL # User wants to do tuning of nFeatures. @@ -548,21 +571,9 @@ setMethod("show", "SelectParams", function(object) { cat("An object of class 'SelectParams'.\n") - IDs <- c("Selection Name", "Ensemble Selection") - index <- na.omit(match(IDs, object@characteristics[, "characteristic"])) - selectText <- object@characteristics[index, "characteristic"] - cat(selectText, ": ", object@characteristics[index, "value"], ".\n", sep = '') - if(selectText == "Ensemble Selection") - cat("Minimum Functions Selected By:", object@minPresence) - - otherInfo <- object@characteristics[-index, ] - if(nrow(otherInfo) > 0) - { - for(rowIndex in 1:nrow(otherInfo)) - { - cat(otherInfo[rowIndex, "characteristic"], ": ", otherInfo[rowIndex, "value"], ".\n", sep = '') - } - } + .showCharacteristics(object@characteristics, c("Selection Name", "Ensemble Selection")) + if(is.list(object@featureRanking)) + cat("Minimum Functions Selected By: ", object@minPresence, ".\n", sep = '') }) ##### TrainParams ##### @@ -588,15 +599,15 @@ setClass("TrainParams", representation( #' @docType class #' @section Constructor: #' \describe{ -#' \item{\code{TrainParams(classifier, balancing = c("downsample", "upsample", "none"), characteristics = DataFrame(), +#' \item{\code{TrainParams(classifier, characteristics = DataFrame(), #' intermediate = character(0), tuneParams = NULL, getFeatures = NULL, ...)}}{ #' Creates a \code{TrainParams} object which stores the function which will do the #' classifier building and parameters that the function will use. #' \describe{ #' \item{\code{classifier}}{A character keyword referring to a registered classifier. See \code{\link{available}} #' for valid keywords.} -#' \item{\code{balancing}}{Default: \code{"downsample"}. A keyword specifying how to handle class imbalance for data sets with categorical outcome. -#' Valid values are \code{"downsample"}, \code{"upsample"} and \code{"none"}.} +#' \item{\code{balancing}}{Not used. Class imbalance is handled by the \code{balancing} setting of \code{\link{ModellingParams}}; +#' specifying it here is an error.} #' \item{\code{characteristics}}{A \code{\link{DataFrame}} describing the #' characteristics of the classifier used. First column must be named \code{"charateristic"} #' and second column must be named \code{"value"}. If using wrapper functions for classifiers @@ -635,12 +646,16 @@ setClassUnion("characterOrFunction", c("character", "function")) #' @rdname TrainParams-class #' @export setMethod("TrainParams", c("characterOrFunction"), - function(classifier, balancing = c("downsample", "upsample", "none"), characteristics = DataFrame(), intermediate = character(0), tuneParams = NULL, getFeatures = NULL, ...) + function(classifier, balancing, characteristics = DataFrame(), intermediate = character(0), tuneParams = NULL, getFeatures = NULL, ...) { + if(!missing(balancing)) + stop("'balancing' is a setting of ModellingParams, not TrainParams. Specify it with ModellingParams(balancing = ...).") extras <- list(...) if(is.character(classifier)) { TPparams <- .classifierKeywordToParams(classifier, tuneParams) + if(is.null(TPparams)) + stop("'", classifier, "' is not a classifier keyword. See available() for valid keywords.") trainParams <- TPparams[[1]] # Get a default params object. if(is.null(getFeatures) && !is.null(trainParams@getFeatures)) getFeatures <- trainParams@getFeatures @@ -650,7 +665,7 @@ setMethod("TrainParams", c("characterOrFunction"), } if(ncol(characteristics) == 0 || !"Classifier Name" %in% characteristics[, "characteristic"]) { - characteristics <- rbind(characteristics, S4Vectors::DataFrame(characteristic = "Classifier Name", value = .ClassifyRenvir[["functionsTable"]][.ClassifyRenvir[["functionsTable"]][, "character"] == attr(classifier, "name"), "name"])) + characteristics <- rbind(characteristics, S4Vectors::DataFrame(characteristic = "Classifier Name", value = .functionDisplayName(classifier, "User-specified Classifier"))) } if(length(extras) == 0) extras <- NULL @@ -664,18 +679,7 @@ setMethod("show", "TrainParams", function(object) { cat("An object of class 'TrainParams'.\n") - index <- na.omit(match("Classifier Name", object@characteristics[, "characteristic"])) - trainText <- object@characteristics[index, "characteristic"] - cat("Classifier Name: ", object@characteristics[index, "value"], ".\n", sep = '') - - otherInfo <- object@characteristics[-index, ] - if(nrow(otherInfo) > 0) - { - for(rowIndex in 1:nrow(otherInfo)) - { - cat(otherInfo[rowIndex, "characteristic"], ": ", otherInfo[rowIndex, "value"], ".\n", sep = '') - } - } + .showCharacteristics(object@characteristics, "Classifier Name") }) ##### PredictParams ##### @@ -712,7 +716,8 @@ setClass("PredictParams", representation( #' function will use. If the training function also makes predictions, this #' must be set to \code{NULL}.} #' \describe{ \item{\code{predictor}}{A character keyword referring to a registered classifier. See \code{\link{available}} -#' for valid keywords.} +#' for valid keywords. Classifiers which train and predict in one function (\code{"kNN"} and \code{"naiveBayes"}) have no +#' separate predictor; for them, use \code{predictParams = NULL} in \code{\link{ModellingParams}}.} #' \item{\code{characteristics}}{A \code{\link{DataFrame}} describing #' the characteristics of the predictor function used. First column must be #' named \code{"charateristic"} and second column must be named \code{"value"}.} @@ -742,8 +747,16 @@ standardGeneric("PredictParams")) setMethod("PredictParams", c("characterOrFunction"), function(predictor, characteristics = DataFrame(), intermediate = character(0), ...) { - if(is.character(predictor)) - predictor <- .classifierKeywordToParams(keyword = predictor, NULL)[[2]]@predictor # Prediction function. + if(is.character(predictor)) + { + keywordParams <- .classifierKeywordToParams(keyword = predictor, NULL) + if(is.null(keywordParams)) + stop("'", predictor, "' is not a classifier keyword. See available() for valid keywords.") + if(is.null(keywordParams[[2]])) + stop("Classifier '", predictor, "' trains and predicts in one function, so it has no separate predictor.\n", + "Use predictParams = NULL in ModellingParams.") + predictor <- keywordParams[[2]]@predictor # Prediction function. + } others <- list(...) if(length(others) == 0) others <- NULL new("PredictParams", predictor = predictor, characteristics = characteristics, @@ -755,23 +768,9 @@ setMethod("show", "PredictParams", function(object) { cat("An object of class 'PredictParams'.\n") - if(ncol(object@characteristics) > 0) - { - index <- na.omit(match("Predictor Name", object@characteristics[, "characteristic"])) - if(length(index) > 0) - cat("Predictor Name: ", object@characteristics[index, "value"], ".\n", sep = '') - - otherInfo <- object@characteristics[-index, ] - if(nrow(otherInfo) > 0) - { - for(rowIndex in 1:nrow(otherInfo)) - { - cat(otherInfo[rowIndex, "characteristic"], ": ", otherInfo[rowIndex, "value"], ".\n", sep = '') - } - } - } - if(is.null(object@predictor)) - cat("Prediction is done by function specified to TrainParams.\n") + if(!is.null(attr(object@predictor, "name"))) + cat("Predictor Function: ", attr(object@predictor, "name"), ".\n", sep = '') + .showCharacteristics(object@characteristics, "Predictor Name") }) setClassUnion("PredictParamsOrNULL", c("PredictParams", "NULL")) diff --git a/R/constants.R b/R/constants.R index 6c12296..bafe651 100644 --- a/R/constants.R +++ b/R/constants.R @@ -12,7 +12,7 @@ "Macro Recall", "higher", "Macro F1", "higher", "Matthews Correlation Coefficient", "higher", - "AUC", "higer", + "AUC", "higher", "C-index", "higher"), ncol = 2, byrow = TRUE, dimnames = list(NULL, c("type", "better")) ) |> as.data.frame() diff --git a/R/crissCrossValidate.R b/R/crissCrossValidate.R index 6b5d598..a243097 100644 --- a/R/crissCrossValidate.R +++ b/R/crissCrossValidate.R @@ -25,9 +25,14 @@ #' @param runTOP Default: \code{FALSE}. If \code{TRUE}, perform the Transferable Omics Prediction (TOP) procedure in a leave-one-dataset-out manner. #' @param verbose Default: 0. A number between 0 and 3 for the amount of progress messages to give. A higher number will produce more messages. #' +#' @details Only the features measured in every data set are used. For \code{trainType = "modelTrain"}, the diagonal of +#' the performance matrix is resubstitution: the model is evaluated on the data set it was trained on. For +#' \code{trainType = "modelTest"}, the diagonal is the cross-validated performance of feature selection and +#' training within that data set. If \code{nFeatures} has several values, the random feature baseline uses the largest one. #' @return A list with elements \code{"real"} for the matrix of pairwise performance metrics using real #' feature selection, \code{"random"} if \code{doRandomFeatures} is \code{TRUE} for metrics of random selection, -#' \code{"top"} if \code{runTOP} is \code{TRUE}, and \code{"params"} for a list of parameters used. +#' \code{"top"} if \code{runTOP} is \code{TRUE}, and \code{"params"} for a list of parameters used. The element +#' \code{"diagonal"} of \code{"params"} states how the diagonal was evaluated. #' #' @author Harry Robertson #' @export @@ -45,28 +50,35 @@ crissCrossValidate <- function(measurements, outcomes, nCores = 1, verbose = 0) { - if(!requireNamespace("TOP", quietly = TRUE)) + if(runTOP && !requireNamespace("TOP", quietly = TRUE)) stop("The package 'TOP' could not be found. Please install it.") trainType <- match.arg(trainType) - extraParams <- NULL - if(length(nFeatures) > 1) # Tune the choice of top number of features. - { - extraParams <- list(tuneCross = list(tuneMode = selectionOptimisation, performanceType = performanceType)) - } - if(!is.list(measurements)) stop("'measurements' is not of type list but is of type", class(measurements)) + if(!is.list(measurements)) stop("'measurements' is not of type list but is of type ", class(measurements)[1]) if(is.null(names(measurements))) stop("Each element of 'measurements' must be named by the name of the data set.") - if(!is.list(outcomes)) stop("'outcomes' is not of type list but is of type", class(outcomes)) + if(!is.list(outcomes)) stop("'outcomes' is not of type list but is of type ", class(outcomes)[1]) + if(length(outcomes) != length(measurements)) stop("'outcomes' must have one element for each data set in 'measurements'.") - isCategorical <- is.character(outcomes[[1]]) && - (length(outcomes[[1]]) == 1 || length(outcomes[[1]]) == nrow(measurements[[1]])) || - is.factor(outcomes[[1]]) + categorical <- mapply(function(measurementsOne, outcomesOne) + { + is.character(outcomesOne) && (length(outcomesOne) == 1 || length(outcomesOne) == nrow(measurementsOne)) || + is.factor(outcomesOne) + }, measurements, outcomes) + if(length(unique(categorical)) > 1) + stop("Some data sets have a categorical outcome and others have a survival outcome. All must be the same kind.") + isCategorical <- categorical[1] # If user didn't specify performanceType, choose "Balanced Accuracy" for classification else "C-index" for survival if(performanceType == "auto") { if(isCategorical) performanceType <- "Balanced Accuracy" else performanceType <- "C-index" } + + extraParams <- NULL + if(length(nFeatures) > 1) # Tune the choice of top number of features. + { + extraParams <- list(tuneCross = list(tuneMode = selectionOptimisation, performanceType = performanceType)) + } # If user left selectionMethod as "auto", pick t-test for categorical or CoxPH for survival if(length(selectionMethod) == 1 && selectionMethod == "auto") { if(isCategorical) selectionMethod <- "t-test" else selectionMethod <- "CoxPH" @@ -76,6 +88,9 @@ crissCrossValidate <- function(measurements, outcomes, if(isCategorical) classifier <- "randomForest" else classifier <- "CoxPH" } + # Models are trained on one data set and applied to the others, so only features measured in every data set are used. + measurements <- .keepSharedFeatures(measurements) + # Keep a copy of the original measurements/outcomes if runTOP is requested if (runTOP) { top_measurements <- measurements @@ -86,7 +101,7 @@ crissCrossValidate <- function(measurements, outcomes, prepareData(measurementsOne, outcomesOne) }, measurements, outcomes, SIMPLIFY = FALSE) - measurements <- lapply(dataCleaned, "[[", 1) + measurements <- .keepSharedFeatures(lapply(dataCleaned, "[[", 1)) outcomes <- lapply(dataCleaned, "[[", 2) if (trainType == "modelTrain") { @@ -107,7 +122,7 @@ crissCrossValidate <- function(measurements, outcomes, # Predict on each dataset performanceAllPairs <- lapply(trainedModels, function(trainedModel) { mapply(function(testData, testOutcomes) { - predictions <- predict(trainedModel, testData[, attr(trainedModel, "featuresForTrain")], outcome = NULL, verbose = verbose) + predictions <- predict(trainedModel, testData[, attr(trainedModel, "featuresForTrain"), drop = FALSE], outcome = NULL, verbose = verbose) if (performanceType == "AUC") { # Must have columns named after each factor level for multi-class AUC @@ -165,9 +180,12 @@ crissCrossValidate <- function(measurements, outcomes, # Build cross-validation parameters crossValParams <- generateCrossValParams(nRepeats, nFolds, nCores, extraParams) - # Evaluate each "trainedModel" on all datasets - performanceAllPairs <- lapply(trainedModels, function(trainedModel) { - mapply(function(measurementsOne, outcomesOne) { + # Evaluate each "trainedModel" on all datasets. On the diagonal, the features chosen in a data set + # would be evaluated on the same samples, so the cross-validation of that data set is used instead. + performanceAllPairs <- mapply(function(trainedModel, trainName) { + mapply(function(measurementsOne, outcomesOne, testName) { + if(testName == trainName) + return(mean(performance(calcCVperformance(trainedModel, performanceType))[[performanceType]])) classifierParams <- .classifierKeywordToParams(classifier, NULL) modellingParams <- ModellingParams( selectParams = SelectParams("previousSelection", @@ -180,8 +198,8 @@ crissCrossValidate <- function(measurements, outcomes, result <- runTests(measurementsOne, outcomesOne, crossValParams, modellingParams) avgPerf <- mean(performance(calcCVperformance(result, performanceType))[[performanceType]]) avgPerf - }, measurements, outcomes, SIMPLIFY = FALSE) - }) + }, measurements, outcomes, names(measurements), SIMPLIFY = FALSE) + }, trainedModels, names(trainedModels), SIMPLIFY = FALSE) realPerformance <- matrix( unlist(performanceAllPairs), @@ -197,16 +215,17 @@ crissCrossValidate <- function(measurements, outcomes, if (doRandomFeatures) { message("Starting random feature selection procedure.") - # For each dataset, pick random features (nFeatures) + # For each dataset, pick random features. If several values of nFeatures are tuned, the largest is used. + nRandom <- min(max(nFeatures), ncol(measurements[[1]])) randomFeatures <- lapply(measurements, function(dataset) { - sample(colnames(dataset), nFeatures) + sample(colnames(dataset), nRandom) }) performanceAllPairs <- lapply(randomFeatures, function(randomFeaturesSet) { mapply(function(testData, testOutcomes) { resultRand <- crossValidate(testData[, randomFeaturesSet, drop = FALSE], testOutcomes, - nFeatures = nFeatures, + nFeatures = nRandom, selectionMethod = "none", classifier = classifier, multiViewMethod = "none", @@ -289,21 +308,38 @@ crissCrossValidate <- function(measurements, outcomes, trainType = trainType, performanceType = performanceType, doRandomFeatures = doRandomFeatures, - runTOP = runTOP) + runTOP = runTOP, + diagonal = if(trainType == "modelTrain") "resubstitution" else "cross-validation") result } +# Subset every data set to the features measured in all of them. +.keepSharedFeatures <- function(measurements) +{ + sharedFeatures <- Reduce(intersect, lapply(measurements, colnames)) + if(length(sharedFeatures) == 0) + stop("The data sets have no features in common.") + nUnshared <- sapply(measurements, ncol) - length(sharedFeatures) + if(any(nUnshared > 0)) + message("Using the ", length(sharedFeatures), " features present in every data set. Features not used per data set: ", + paste(names(measurements), nUnshared, sep = " ", collapse = ", "), ".") + lapply(measurements, function(measurementsOne) measurementsOne[, sharedFeatures, drop = FALSE]) +} + #' A function to plot the output of the crissCrossValidate function. #' #' This function generates a heatmap of the cross-validation results from -#' \code{\link{crissCrossValidate}}. By default, it hides the "resubstitution" diagonal -#' (where the training == test set) unless \code{showResubMetric = TRUE}. +#' \code{\link{crissCrossValidate}}. For \code{trainType = "modelTrain"}, the diagonal +#' (where the training set is the test set) is resubstitution and is hidden unless \code{showResubMetric = TRUE}. +#' For \code{trainType = "modelTest"}, the diagonal is cross-validation within each data set; it is shown and +#' outlined in black. #' #' @param crissCrossResult The output of the \code{\link{crissCrossValidate}} function. #' @param includeValues Logical. If \code{TRUE}, numeric values are printed on each tile. -#' @param showResubMetric Logical. If \code{FALSE}, the diagonal (resubstitution) cells +#' @param showResubMetric Logical. If \code{FALSE}, the diagonal (resubstitution) cells of a \code{"modelTrain"} result #' are set to \code{NA} and appear grayed-out or blank. Defaults to \code{FALSE}. +#' @return A \code{ggplot} object, or a combined plot of two heatmaps if random features were evaluated. #' #' @import ggplot2 #' @import reshape2 @@ -313,22 +349,24 @@ crissCrossPlot <- function(crissCrossResult, includeValues = FALSE, showResubMetric = FALSE) { - # We'll attach for convenience - attach(crissCrossResult) - on.exit(detach(crissCrossResult), add = TRUE) - + params <- crissCrossResult[["params"]] scalebar_title <- params$performanceType + # Results created before the diagonal was recorded have a resubstitution diagonal. + isResubstitution <- !identical(params$diagonal, "cross-validation") + diagonalText <- if(isResubstitution) "resubstitution" else "cross-validation within the data set" # Helper function: turn matrix into heatmap plotMatrix <- function(mat, xlab_text, ylab_text) { # Convert to matrix if needed mat <- as.matrix(mat) - if (!showResubMetric) { + if (isResubstitution && !showResubMetric) { diag(mat) <- NA } melted_df <- reshape2::melt(mat, na.rm = FALSE, value.name = "value") + diagonalIndices <- cbind(seq_len(min(dim(mat))), seq_len(min(dim(mat)))) + diagonal_df <- data.frame(Var1 = rownames(mat)[diagonalIndices[, 1]], Var2 = colnames(mat)[diagonalIndices[, 2]]) gg <- ggplot(melted_df, aes(x = Var1, y = Var2, fill = value)) + geom_tile(color = "white") + @@ -351,6 +389,12 @@ crissCrossPlot <- function(crissCrossResult, ) + coord_fixed() + if (!isResubstitution || showResubMetric) { + gg <- gg + geom_tile(data = diagonal_df, aes(x = Var1, y = Var2), inherit.aes = FALSE, + fill = NA, colour = "black", linewidth = 0.5) + + labs(caption = paste("Diagonal:", diagonalText)) + } + if (includeValues) { gg <- gg + geom_text(aes(label = round(value, 3)), color = "black", size = 3, na.rm = TRUE) } @@ -359,20 +403,20 @@ crissCrossPlot <- function(crissCrossResult, if (params$trainType == "modelTrain") { - mainMatrix <- real - if ("top" %in% names(crissCrossResult)) { # If runTOP is TRUE and we have result$top, append it as an extra row - mainMatrix <- rbind(mainMatrix, crissCrossResult$top) + mainMatrix <- crissCrossResult[["real"]] + if (!is.null(crissCrossResult[["top"]])) { # If runTOP is TRUE and we have result$top, append it as an extra row + mainMatrix <- rbind(mainMatrix, crissCrossResult[["top"]]) } heatmapObj <- plotMatrix(mainMatrix, "Training Dataset", "Testing Dataset") } else if (params$trainType == "modelTest") { # Real is "Features Extracted" vs. "Cross-validate" - mainMatrix <- real + mainMatrix <- crissCrossResult[["real"]] heatmapObj1 <- plotMatrix(mainMatrix, "Features Extracted", "Dataset Tested") - if (params$doRandomFeatures == TRUE && exists("random")) { - heatmapObj2 <- plotMatrix(random, "Random Features", "Dataset Tested") + if (isTRUE(params$doRandomFeatures) && !is.null(crissCrossResult[["random"]])) { + heatmapObj2 <- plotMatrix(crissCrossResult[["random"]], "Random Features", "Dataset Tested") heatmapObj <- ggarrange( heatmapObj1, heatmapObj2, @@ -386,5 +430,5 @@ crissCrossPlot <- function(crissCrossResult, } } - print(heatmapObj) + heatmapObj } diff --git a/R/data.R b/R/data.R index 834bf4a..8385774 100644 --- a/R/data.R +++ b/R/data.R @@ -38,13 +38,16 @@ NULL #' METABRIC Clinical Data #' -#' 470 patients with eight features. +#' 470 patients with six clinical features and their relapse-free survival. #' #' @name METABRICclinical #' @aliases clinical #' @docType data -#' @format \code{clinical} A \code{\link{DataFrame}} containing clinical data. +#' @format \code{clinical} A \code{\link{DataFrame}} with a row for each patient. The six features are +#' \code{Breast.Tumour.Laterality}, \code{ER.Status}, \code{Inferred.Menopausal.State}, \code{Lymph.Nodes.Positive}, +#' \code{Grade} and \code{Size}. The outcome is relapse-free survival: \code{timeRFS} is the time and +#' \code{eventRFS} is the event indicator (1 for relapse). Use them as the outcome, not as features. #' @source Dynamics of Breast Cancer Relapse Reveal Late-recurring ER-positive Genomic Subgroups, \emph{Nature}, -#' 2019. Webpage: \url{https://www.nature.com/articles/s43018-020-0026-6} +#' 2019. Webpage: \url{https://www.nature.com/articles/s41586-019-1007-8} #' @keywords datasets NULL diff --git a/R/interfaceRandomForest.R b/R/interfaceRandomForest.R index 2979741..3a5f056 100644 --- a/R/interfaceRandomForest.R +++ b/R/interfaceRandomForest.R @@ -17,6 +17,8 @@ attr(randomForestTrainInterface, "name") <- "randomForestTrainInterface" # forest is of class ranger randomForestPredictInterface <- function(forest, measurementsTest, ..., returnType = c("both", "class", "score"), verbose = 3) { + if(!requireNamespace("ranger", quietly = TRUE)) + stop("The package 'ranger' could not be found. Please install it.") returnType <- match.arg(returnType) classes <- forest$forest$levels if(verbose == 3) diff --git a/R/precisionPathways.R b/R/precisionPathways.R index dde70c7..8bd224e 100644 --- a/R/precisionPathways.R +++ b/R/precisionPathways.R @@ -7,7 +7,7 @@ #' #' @param measurements Either a \code{\link{MultiAssayExperiment}} or a list of the basic tabular objects containing the data. #' @param class If a \code{\link{MultiAssayExperiment}}, a column name in \code{colData(measurements)} with the classes. If \code{measurements} is a \code{list} of tabular data, may also be -#' a vector of classes. +#' a vector of classes. Optional for \code{precisionPathwaysPredict}: the classes of the new samples are only needed to calculate accuracies. #' @param useFeatures Default: \code{NULL} (i.e. use all provided features). A named list of features to use. Otherwise, the input data is a single table and this can just be a vector of feature names. #' For any assays not in the named list, all of their features are used. \code{"clinical"} is also a valid assay name and refers to the clinical data table. #' This allows for the avoidance of variables such spike-in RNAs, sample IDs, sample acquisition dates, etc. which are not relevant for outcome prediction. @@ -16,9 +16,10 @@ #' @param topNvariance Default: NULL. An integer number of most variable features per assay to subset to. #' Assays with less features won't be reduced in size. #' @param fixedAssays A character vector of assay names specifying any assays which must be at the -#' beginning of the pathway. +#' beginning of the pathway, in the order given. \code{NULL} allows any assay to be first. #' @param confidenceCutoff The minimum confidence of predictions for a sample to be predicted by a particular assay. If a sample #' was predicted to belong to a particular class a proportion \eqn{p} times, then the confidence is \eqn{2 \times |p - 0.5|}. +#' A sample whose confidence is equal to or greater than \code{confidenceCutoff} is predicted at that tier. #' @param minAssaySamples An integer specifying the minimum number of samples a tier may have. If a subsequent tier #' would have less than this number of samples, the samples are incorporated into the current tier. #' @param mode Default: \code{"stability"}. Either \code{"stability"} or \code{"combinatorial"}. If \code{"stability"}, then pathways grow by @@ -32,6 +33,12 @@ #' @param nRepeats A numeric specifying the the number of repeats or permutations to use for cross-validation. #' @param nCores A numeric specifying the number of cores used if the user wants to use parallelisation. #' @param pathways A set of pathways created by \code{precisionPathwaysTrain} which is an object of class \code{PrecisionPathways} to be used for predicting on a new data set. +#' @details \code{precisionPathwaysTrain} cross-validates a classifier for each assay. The confidence of a sample for an assay is +#' calculated from its cross-validated predictions. \code{precisionPathwaysPredict} predicts each new sample with every model fitted +#' during that cross-validation (one model per fold of every repeat), and the confidence of a new sample for an assay is calculated from +#' those predictions. Each new sample is predicted independently of the classes of the new samples; the stored models are used as they are. +#' The pathways are then formed in the same way as in training. Each pathway records the assays in the order in which samples passed through +#' them, in its element named \code{"assays"}. #' @rdname precisionPathways #' @aliases precisionPathwaysTrain precisionPathwaysPredict #' @return An object of class \code{PrecisionPathways} which is basically a named list that other plotting and @@ -45,7 +52,7 @@ setGeneric("precisionPathwaysTrain", function(measurements, class, ...) #' @rdname precisionPathways #' @export -setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", +setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", function(measurements, class, useFeatures = NULL, maxMissingProp = 0.0, topNvariance = NULL, fixedAssays = "clinical", confidenceCutoff = 0.8, minAssaySamples = 10, mode = c("stability", "combinatorial"), nFeatures = 20, selectionMethod = setNames(c("none", rep("t-test", length(measurements))), c("clinical", names(measurements))), @@ -57,19 +64,28 @@ setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", # One of the tables must be named "clinical". if (!any(names(measurements) == "clinical")) stop("One of the tables must be named \"clinical\".") + clinicalColumns <- colnames(measurements[["clinical"]]) + } else { + clinicalColumns <- colnames(MultiAssayExperiment::colData(measurements)) } - if(is.null(useFeatures) && fixedAssays == "clinical") + # A list input already has a table named "clinical", so the default vectors name it twice. + selectionMethod <- selectionMethod[!duplicated(names(selectionMethod))] + classifier <- classifier[!duplicated(names(classifier))] + + if(is.null(useFeatures) || !"clinical" %in% names(useFeatures)) { - warning("No 'useFeatures' named list element for clinical data is specified. Clinical data often has\n", - "lots of uninformative variables. Please consider specifying useful features.") - useFeatures <- list(clinical = colnames(MultiAssayExperiment::colData(measurements))) + warning("No 'useFeatures' named list element for clinical data is specified. Clinical data often has\n", + "lots of uninformative variables. Please consider specifying useful features.") + if(is.character(class)) clinicalColumns <- setdiff(clinicalColumns, class) + if(is.null(useFeatures)) useFeatures <- list() + useFeatures[["clinical"]] <- clinicalColumns } mode <- match.arg(mode) - - prepArgs <- list(measurements, outcomeColumns = class, useFeatures = useFeatures, + + prepArgs <- list(measurements, class, useFeatures = useFeatures, maxMissingProp = maxMissingProp, topNvariance = topNvariance) measurementsAndClass <- do.call(prepareData, prepArgs) - + .precisionPathwaysTrain(measurementsAndClass[["measurements"]], measurementsAndClass[["outcome"]], useFeatures = useFeatures, fixedAssays = fixedAssays, confidenceCutoff = confidenceCutoff, minAssaySamples = minAssaySamples, mode = mode, nFeatures = nFeatures, @@ -88,8 +104,13 @@ setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", # Step 1: Determine all valid permutations of assays, taking into account the # assays to be used and which assays, if any, must be included. assayIDs <- unique(S4Vectors::mcols(measurements)[["assay"]]) - assaysPermutations <- .permutations(assayIDs, fixed = data.frame(seq_along(fixedAssays), fixedAssays)) - assaysPermutations <- asplit(assaysPermutations, 2) + if(!all(fixedAssays %in% assayIDs)) + stop("'fixedAssays' ", paste(setdiff(fixedAssays, assayIDs), collapse = ", "), " not found in the data. Assays are: ", + paste(assayIDs, collapse = ", "), '.') + fixedPositions <- NULL + if(length(fixedAssays) > 0) fixedPositions <- data.frame(seq_along(fixedAssays), fixedAssays) + assaysPermutations <- .permutations(assayIDs, fixed = fixedPositions) + assaysPermutations <- lapply(asplit(assaysPermutations, 2), as.vector) if(mode == "combinatorial") { subsetPermutations <- list() @@ -111,83 +132,138 @@ setMethod("precisionPathwaysTrain", "MultiAssayExperimentOrList", assaysPermutations <- c(assaysPermutations, assayIDs) # Each assay on its own. } - permutationIDs <- lapply(assaysPermutations, function(permutation) paste(permutation, collapse = '-')) - # Step 2: Build a classifier for each assay using all of the samples. modelsList <- crossValidate(measurements, class, nFeatures, selectionMethod, classifier = classifier, nFolds = nFolds, nRepeats = nRepeats, nCores = nCores) + if(is(modelsList, "ClassifyResult")) modelsList <- list(modelsList) modelsList <- lapply(modelsList, calcCVperformance, "Sample Accuracy") # Add sample accuracy, which can be subset later. + names(modelsList) <- .pathwaysAssayNames(modelsList, assayIDs) + modelsList <- lapply(modelsList, .pathwaysTrimModels) # Step 3: Loop over each pathway and each assay in order to determine which samples are used at that level # and which are passed onwards. - precisionPathways <- lapply(assaysPermutations, function(permutation) - { - assaysProcessed <- character() - samplesUsed <- character() - individualsTableAll <- S4Vectors::DataFrame() - tierTableAll <- S4Vectors::DataFrame() - breakEarly = FALSE - for(assay in permutation) - { - # Step 3a: Identify all samples which are consistently predicted. - modelIndex <- match(assay, assayIDs) - allPredictions <- predictions(modelsList[[modelIndex]]) - allSampleIDs <- sampleNames(modelsList[[modelIndex]]) - predictionsSamplesCounts <- table(allPredictions[, "sample"], allPredictions[, "class"]) - confidences <- 2 * abs(predictionsSamplesCounts[, 1] / rowSums(predictionsSamplesCounts) - 0.5) - sampleIDsUse <- names(confidences)[confidences > confidenceCutoff] - sampleIDsUse <- setdiff(sampleIDsUse, samplesUsed) - - # Check if too few samples left for next round. Include them in this round, if so. - remainingIDs <- setdiff(allSampleIDs, c(samplesUsed, sampleIDsUse)) - if(length(remainingIDs) < minAssaySamples || assay == permutation[length(permutation)]) - { - sampleIDsUse <- c(sampleIDsUse, remainingIDs) - breakEarly = TRUE - } - if(length(sampleIDsUse) == 0) next # No samples have sufficient confidence at this tier. - - predictionsSamplesCounts <- predictionsSamplesCounts[sampleIDsUse, , drop = FALSE] - - # Step 3b: Individuals predictions and sample-wise accuracy, tier-wise error. - maxVotes <- apply(predictionsSamplesCounts, 1, function(sample) which.max(sample)) - predictedClasses <- factor(colnames(predictionsSamplesCounts)[maxVotes], - levels = colnames(predictionsSamplesCounts)) - individualsTable <- S4Vectors::DataFrame(Tier = assay, - `Sample ID` = sampleIDsUse, - `Predicted` = predictedClasses, - `Accuracy` = performance(modelsList[[modelIndex]])[["Sample Accuracy"]][sampleIDsUse], - check.names = FALSE) - knownClasses <- actualOutcome(modelsList[[modelIndex]])[match(sampleIDsUse, allSampleIDs)] - balancedAccuracy <- calcExternalPerformance(knownClasses, predictedClasses) - tierTable <- S4Vectors::DataFrame(Tier = assay, - `Balanced Accuracy` = balancedAccuracy, check.names = FALSE) - - assaysProcessed <- c(assaysProcessed, assay) - individualsTableAll <- rbind(individualsTableAll, individualsTable) - tierTableAll <- rbind(tierTableAll, tierTable) - samplesUsed <- c(samplesUsed, sampleIDsUse) - - if(breakEarly == TRUE) break - } - pathwayString <- paste(assaysProcessed, collapse = '-') - individualsTableAll[, "Tier"] <- factor(individualsTableAll[, "Tier"], levels = permutation) - list(pathway = pathwayString, - individuals = individualsTableAll, tiers = tierTableAll) - }) + predictionsSamplesCounts <- lapply(modelsList, function(result) + { + allPredictions <- predictions(result) + table(allPredictions[, "sample"], allPredictions[, "class"]) + }) + sampleAccuracies <- lapply(modelsList, function(result) performance(result)[["Sample Accuracy"]]) + allSampleIDs <- sampleNames(modelsList[[1]]) + knownClasses <- actualOutcome(modelsList[[1]]) + precisionPathways <- lapply(assaysPermutations, .pathwayTiers, predictionsSamplesCounts, sampleAccuracies, + allSampleIDs, knownClasses, confidenceCutoff, minAssaySamples) names(precisionPathways) <- sapply(precisionPathways, "[[", "pathway") - attr(modelsList, "crossValParams") <- generateCrossValParams(nRepeats, nFolds, nCores, NULL) + attr(modelsList, "crossValParams") <- generateCrossValParams(nRepeats, nFolds, nCores, NULL) precisionPathways <- precisionPathways[unique(names(precisionPathways))] # In case early termination causes duplicates. result <- list(models = modelsList, assaysPermutations = assaysPermutations, parameters = list(confidenceCutoff = confidenceCutoff, minAssaySamples = minAssaySamples, classifier = classifier), useFeatures = useFeatures, pathways = precisionPathways) class(result) <- "PrecisionPathways" - + result } +# The assay name of each cross-validation result, so that models are always looked up by assay and never by position. +.pathwaysAssayNames <- function(modelsList, assayIDs) +{ + assayNames <- sapply(modelsList, function(result) + { + characteristics <- result@characteristics + assayName <- characteristics[characteristics[, "characteristic"] == "Assay Name", "value"] + if(length(assayName) == 0) NA_character_ else assayName + }) + if(anyNA(assayNames) && length(modelsList) == length(assayIDs)) # Results lacking an assay name are in data order. + assayNames[is.na(assayNames)] <- assayIDs[is.na(assayNames)] + assayNames +} + +# Prediction of new samples uses only the fitted models. The extra forests that some classifiers +# store for calculating variable importance are not needed. +.pathwaysTrimModels <- function(result) +{ + result@models <- lapply(result@models, function(model) + { + attr(model, "forImportance") <- NULL + model + }) + result +} + +# Confidence of each sample, from a table of counts of predicted classes (samples in rows, classes in columns). +# A small tolerance ensures that a confidence equal to the cut-off is not missed because of floating point error. +.pathwaysConfident <- function(predictionsSamplesCounts, confidenceCutoff) +{ + confidences <- 2 * abs(predictionsSamplesCounts[, 1] / rowSums(predictionsSamplesCounts) - 0.5) + names(confidences)[confidences >= confidenceCutoff - sqrt(.Machine$double.eps)] +} + +# Balanced accuracy, with predicted and known classes using the same factor levels in the same order. +.pathwaysBalancedAccuracy <- function(knownClasses, predictedClasses) +{ + if(is.null(knownClasses) || all(is.na(knownClasses))) return(NA_real_) + classLevels <- union(levels(predictedClasses), levels(knownClasses)) + calcExternalPerformance(factor(as.character(knownClasses), levels = classLevels), + factor(as.character(predictedClasses), levels = classLevels)) +} + +# Pass samples along one sequence of assays. At each tier, the samples which have not yet been predicted and +# are predicted with sufficient confidence are predicted by that tier. Used for training and for prediction. +.pathwayTiers <- function(permutation, predictionsSamplesCounts, sampleAccuracies, allSampleIDs, knownClasses, + confidenceCutoff, minAssaySamples) +{ + assaysProcessed <- character() + assaysTested <- character() + samplesUsed <- character() + individualsTableAll <- S4Vectors::DataFrame() + tierTableAll <- S4Vectors::DataFrame() + breakEarly = FALSE + for(assay in permutation) + { + # Identify all samples which are consistently predicted. + predictionsSamplesCountsUse <- predictionsSamplesCounts[[assay]] + sampleIDsUse <- .pathwaysConfident(predictionsSamplesCountsUse, confidenceCutoff) + sampleIDsUse <- setdiff(sampleIDsUse, samplesUsed) + assaysTested <- c(assaysTested, assay) # Every sample not yet predicted has this assay done. + + # Check if too few samples left for next round. Include them in this round, if so. + remainingIDs <- setdiff(allSampleIDs, c(samplesUsed, sampleIDsUse)) + if(length(remainingIDs) < minAssaySamples || assay == permutation[length(permutation)]) + { + sampleIDsUse <- c(sampleIDsUse, remainingIDs) + breakEarly = TRUE + } + if(length(sampleIDsUse) == 0) next # No samples have sufficient confidence at this tier. + + predictionsSamplesCountsUse <- predictionsSamplesCountsUse[sampleIDsUse, , drop = FALSE] + + # Individuals predictions and sample-wise accuracy, tier-wise error. + maxVotes <- apply(predictionsSamplesCountsUse, 1, function(sample) which.max(sample)) + predictedClasses <- factor(colnames(predictionsSamplesCountsUse)[maxVotes], + levels = colnames(predictionsSamplesCountsUse)) + individualsTable <- S4Vectors::DataFrame(Tier = assay, + `Sample ID` = sampleIDsUse, + `Predicted` = predictedClasses, + `Accuracy` = sampleAccuracies[[assay]][sampleIDsUse], + check.names = FALSE) + balancedAccuracy <- .pathwaysBalancedAccuracy(knownClasses[match(sampleIDsUse, allSampleIDs)], predictedClasses) + tierTable <- S4Vectors::DataFrame(Tier = assay, + `Balanced Accuracy` = balancedAccuracy, check.names = FALSE) + + assaysProcessed <- c(assaysProcessed, assay) + individualsTableAll <- rbind(individualsTableAll, individualsTable) + tierTableAll <- rbind(tierTableAll, tierTable) + samplesUsed <- c(samplesUsed, sampleIDsUse) + + if(breakEarly == TRUE) break + } + pathwayString <- paste(assaysProcessed, collapse = '-') + individualsTableAll[, "Tier"] <- factor(individualsTableAll[, "Tier"], levels = permutation) + list(pathway = pathwayString, assays = assaysProcessed, assaysTested = assaysTested, + individuals = individualsTableAll, tiers = tierTableAll) +} + # A nice print method to avoid flooding the screen with lots of tables # when result is shown in console. #' @export @@ -196,6 +272,8 @@ print.PrecisionPathways <- function(x, ...) cat("An object of class 'PrecisionPathways'.\n") cat("Pathways:\n") cat(paste(names(x[["pathways"]]), collapse = '\n')) + cat("\n") + invisible(x) } #' @usage NULL @@ -204,9 +282,10 @@ setGeneric("precisionPathwaysPredict", function(pathways, measurements, class, . #' @rdname precisionPathways #' @export -setMethod("precisionPathwaysPredict", c("PrecisionPathways", "MultiAssayExperimentOrList"), +setMethod("precisionPathwaysPredict", c("PrecisionPathways", "MultiAssayExperimentOrList"), function(pathways, measurements, class) { + if(missing(class)) class <- NULL # Classes of new samples are only used for calculating accuracy. if(is.list(measurements)) # Ensure plain list has clinical data. { # One of the tables must be named "clinical". @@ -214,138 +293,124 @@ setMethod("precisionPathwaysPredict", c("PrecisionPathways", "MultiAssayExperime stop("One of the tables must be named \"clinical\".") } - prepArgs <- list(measurements, outcomeColumns = class, useFeatures = pathways[["useFeatures"]]) - measurementsAndClass <- do.call(prepareData, prepArgs) - + if(is.null(class) && is(measurements, "MultiAssayExperiment")) + { # prepareData needs an outcome column. Sample identifiers stand in for it and are then discarded. + class <- ".pathwaysSampleID" + MultiAssayExperiment::colData(measurements)[, class] <- rownames(MultiAssayExperiment::colData(measurements)) + measurementsAndClass <- prepareData(measurements, class, useFeatures = pathways[["useFeatures"]]) + measurementsAndClass["outcome"] <- list(NULL) + } else { + prepArgs <- list(measurements, class, useFeatures = pathways[["useFeatures"]]) + measurementsAndClass <- do.call(prepareData, prepArgs) + } + .precisionPathwaysPredict(pathways, measurementsAndClass[["measurements"]], measurementsAndClass[["outcome"]]) }) .precisionPathwaysPredict <- function(pathways, measurements, class) { - # Step 1: Extract all of previously fitted models and permutations. Create new predictions per assay. + # Step 1: Extract all of previously fitted models and pathways. Predict the new samples with every fitted model of each assay. modelsList <- pathways[["models"]] - assayIDs <- lapply(pathways[["models"]], function(model) model@characteristics[model@characteristics[, 1] == "Assay Name", 2]) - assaysPermutations <- strsplit(names(pathways[["pathways"]]), '-') + names(modelsList) <- .pathwaysAssayNames(modelsList, names(modelsList)) + assaysPermutations <- lapply(pathways[["pathways"]], function(pathway) + { + if(!is.null(pathway[["assays"]])) pathway[["assays"]] else as.character(pathway[["tiers"]][, "Tier"]) + }) confidenceCutoff <- pathways[["parameters"]][["confidenceCutoff"]] minAssaySamples <- pathways[["parameters"]][["minAssaySamples"]] classifiers <- pathways[["parameters"]][["classifier"]] allSampleIDs <- rownames(measurements) - - resultsNew <- mapply(function(fittedModel, classifierID) + + resultsNew <- mapply(function(fittedModel, assayID) { - trainParams <- TrainParams("previousTrained", classifyResult = fittedModel, intermediate = ".iteration") - predictPrams <- PredictParams(classifierID) - modelParams <- ModellingParams(selectParams = NULL, trainParams = trainParams, predictParams = predictPrams) - resultNew <- runTests(measurements, class, attr(modelsList, "crossValParams"), modelParams, verbose = 0) - resultNew <- calcCVperformance(resultNew, "Sample Accuracy") - predictsNew <- predictions(resultNew) + if(!is.null(fittedModel@modellingParams)) + predictParams <- fittedModel@modellingParams@predictParams + else # Result created without its modelling parameters. + predictParams <- PredictParams(classifiers[[assayID]]) + if(is.null(predictParams)) + stop("The classifier of assay ", assayID, " trains and predicts in one step, so its models cannot predict new samples.") + classLevels <- levels(actualOutcome(fittedModel)) + predictedClasses <- lapply(models(fittedModel), function(model) + { + measurementsUse <- measurements[, attr(model, "featuresForTrain"), drop = FALSE] + paramList <- c(list(model, measurementsUse), predictParams@otherParams, verbose = 0) + predicted <- do.call(predictParams@predictor, paramList) + if(!is.factor(predicted)) predicted <- predicted[, "class"] + as.character(predicted) + }) + predictsNew <- S4Vectors::DataFrame(sample = rep(allSampleIDs, length(predictedClasses)), + class = factor(unlist(predictedClasses), levels = classLevels)) predictionsSamplesCounts <- table(predictsNew[, "sample"], predictsNew[, "class"]) - sampleAccuracies <- performance(resultNew)[["Sample Accuracy"]] - confidences <- 2 * abs(predictionsSamplesCounts[, 1] / rowSums(predictionsSamplesCounts) - 0.5) - list(predictionsSamplesCounts, confidences, sampleAccuracies) - }, modelsList, classifiers, SIMPLIFY = FALSE) + sampleAccuracies <- setNames(rep(NA_real_, nrow(predictionsSamplesCounts)), rownames(predictionsSamplesCounts)) + if(!is.null(class)) + { + knownClasses <- as.character(class)[match(rownames(predictionsSamplesCounts), allSampleIDs)] + correctColumn <- match(knownClasses, colnames(predictionsSamplesCounts)) + correctCounts <- predictionsSamplesCounts[cbind(seq_len(nrow(predictionsSamplesCounts)), correctColumn)] + correctCounts[is.na(correctCounts)] <- 0 # Known class never predicted by the models. + sampleAccuracies[] <- correctCounts / rowSums(predictionsSamplesCounts) + } + list(predictionsSamplesCounts, sampleAccuracies) + }, modelsList, names(modelsList), SIMPLIFY = FALSE) predictionsSamplesCounts <- lapply(resultsNew, "[[", 1) - confidences <- lapply(resultsNew, "[[", 2) - sampleAccuracies <- lapply(resultsNew, "[[", 3) - names(predictionsSamplesCounts) <- names(confidences) <- names(sampleAccuracies) <- assayIDs + sampleAccuracies <- lapply(resultsNew, "[[", 2) # Step 2: Loop over each pathway and each assay in order to determine which samples are used at that level # and which are passed onwards. - precisionPathways <- lapply(assaysPermutations, function(permutation) - { - assaysProcessed <- character() - samplesUsed <- character() - individualsTableAll <- S4Vectors::DataFrame() - tierTableAll <- S4Vectors::DataFrame() - breakEarly = FALSE - for(assay in permutation) - { - # Step 2a: Identify new samples which are consistently predicted. - confidencesUse <- confidences[[assay]] - predictionsSamplesCountsUse <- predictionsSamplesCounts[[assay]] - sampleIDsUse <- names(confidencesUse)[confidencesUse > confidenceCutoff] - sampleIDsUse <- setdiff(sampleIDsUse, samplesUsed) - - # Check if too few samples left for next round. Include them in this round, if so. - remainingIDs <- setdiff(allSampleIDs, c(samplesUsed, sampleIDsUse)) - if(length(remainingIDs) < minAssaySamples || assay == permutation[length(permutation)]) - { - sampleIDsUse <- c(sampleIDsUse, remainingIDs) - breakEarly = TRUE - } - if(length(sampleIDsUse) == 0) next # No samples have sufficient confidence at this tier. - - predictionsSamplesCountsUse <- predictionsSamplesCountsUse[sampleIDsUse, , drop = FALSE] - - # Step 2b: Individuals predictions and sample-wise accuracy, tier-wise error. - maxVotes <- apply(predictionsSamplesCountsUse, 1, function(sample) which.max(sample)) - predictedClasses <- factor(colnames(predictionsSamplesCountsUse)[maxVotes], - levels = colnames(predictionsSamplesCountsUse)) - individualsTable <- S4Vectors::DataFrame(Tier = assay, - `Sample ID` = sampleIDsUse, - `Predicted` = predictedClasses, - `Accuracy` = sampleAccuracies[[assay]][sampleIDsUse], - check.names = FALSE) - knownClasses <- class[match(sampleIDsUse, allSampleIDs)] - balancedAccuracy <- calcExternalPerformance(knownClasses, predictedClasses) - tierTable <- S4Vectors::DataFrame(Tier = assay, - `Balanced Accuracy` = balancedAccuracy, check.names = FALSE) - - assaysProcessed <- c(assaysProcessed, assay) - individualsTableAll <- rbind(individualsTableAll, individualsTable) - tierTableAll <- rbind(tierTableAll, tierTable) - samplesUsed <- c(samplesUsed, sampleIDsUse) - - if(breakEarly == TRUE) break - } - pathwayString <- paste(assaysProcessed, collapse = '-') - individualsTableAll[, "Tier"] <- factor(individualsTableAll[, "Tier"], levels = permutation) - list(pathway = pathwayString, - individuals = individualsTableAll, tiers = tierTableAll) - }) + precisionPathways <- lapply(assaysPermutations, .pathwayTiers, predictionsSamplesCounts, sampleAccuracies, + allSampleIDs, class, confidenceCutoff, minAssaySamples) names(precisionPathways) <- sapply(precisionPathways, "[[", "pathway") precisionPathways <- precisionPathways[unique(names(precisionPathways))] - result <- list(models = modelsList, assaysPermutations = assaysPermutations, - parameters = list(confidenceCutoff = confidenceCutoff, minAssaySamples = minAssaySamples), - pathways = precisionPathways, testSampleInfo = DataFrame(sampleID = rownames(measurements), class = class)) + if(is.null(class)) class <- factor(rep(NA, length(allSampleIDs)), levels = levels(actualOutcome(modelsList[[1]]))) + result <- list(models = modelsList, assaysPermutations = unname(assaysPermutations), + parameters = list(confidenceCutoff = confidenceCutoff, minAssaySamples = minAssaySamples, classifier = classifiers), + pathways = precisionPathways, testSampleInfo = DataFrame(sampleID = allSampleIDs, class = class)) class(result) <- "PrecisionPathways" - + result } # Calculate accuracy and costs of each pathway. #' Various Functions for Evaluating Precision Pathways -#' +#' #' These functions tabulate or plot various aspects of precision pathways, such as accuracies and costs. -#' +#' +#' The cost of a pathway is the total cost of the assays done on all of the samples. A sample predicted at a particular tier +#' had every assay of the pathway up to and including that tier done, so it is charged the cost of each of those assays. +#' For example, a sample predicted by the second assay of the pathway clinical-RNA costs the price of clinical data plus the price of RNA. +#' #' @param precisionPathways A pathway of class \code{PrecisionPathways}. -#' @param costs A named vector of assays with the cost of each one. +#' @param costs A named vector of assays with the cost of each one. An assay that is not named costs nothing. #' @rdname precisionPathwaysEvaluations #' @export calcCostsAndPerformance <- function(precisionPathways, costs = NULL) { if(is.null(costs)) stop("'costs' of each assay must be specified.") - - pathwayIDs <- names(precisionPathways[["pathways"]]) + if("testSampleInfo" %in% names(precisionPathways)) { knownClasses <- precisionPathways[["testSampleInfo"]][, "class"] allNames <- precisionPathways[["testSampleInfo"]][, "sampleID"] } else{ knownClasses <- actualOutcome(precisionPathways$models[[1]]) - allNames <- sampleNames(precisionPathways$models[[1]]) + allNames <- sampleNames(precisionPathways$models[[1]]) } accuraciesCosts <- do.call(rbind, lapply(precisionPathways[["pathways"]], function(pathway) { predictions <- pathway[["individuals"]][, "Predicted"] knownClasses <- knownClasses[match(pathway[["individuals"]][, "Sample ID"], allNames)] - balancedAccuracy <- calcExternalPerformance(knownClasses, predictions) - costTotal <- sum(costs[na.omit(match(pathway[["individuals"]][, "Tier"], names(costs)))]) - if(is.na(costTotal)) # Only clinical used. - costTotal <- 0 - + balancedAccuracy <- .pathwaysBalancedAccuracy(knownClasses, predictions) + + # Each sample is charged for every assay it passed through, up to and including the tier that predicted it. + assaysTested <- pathway[["assaysTested"]] + if(is.null(assaysTested)) assaysTested <- as.character(pathway[["tiers"]][, "Tier"]) + assaysCosts <- unname(costs[match(assaysTested, names(costs))]) + assaysCosts[is.na(assaysCosts)] <- 0 + cumulativeCosts <- setNames(cumsum(assaysCosts), assaysTested) + costTotal <- sum(cumulativeCosts[as.character(pathway[["individuals"]][, "Tier"])]) + data.frame(accuracy = round(balancedAccuracy, 2), cost = costTotal) })) @@ -356,8 +421,8 @@ calcCostsAndPerformance <- function(precisionPathways, costs = NULL) # Print a summary table, including accuracy and costs. #' @param object A set of pathways of class \code{PrecisionPathways}. -#' @param weights A numeric vector of length two specifying how to weight the predictive accuracy -#' and the cost during ranking. Must sum to 1. +#' @param weights A numeric vector of length two named \code{"accuracy"} and \code{"cost"} specifying how to weight the predictive accuracy +#' and the cost during ranking. Must sum to 1. An unnamed vector is taken to be in the order accuracy, cost. #' @param ... Not used but just following the S3 requirement of the generic template. #' @rdname precisionPathwaysEvaluations #' @export @@ -365,19 +430,22 @@ summary.PrecisionPathways <- function(object, weights = c(accuracy = 0.5, cost = { if(!"performance" %in% names(object)) stop("Nothing to summarise. Please run function 'calcCostsAndPerformance' first.") - hasCost <- "cost" %in% names(object[["performance"]]) + if(is.null(names(weights))) names(weights) <- c("accuracy", "cost")[seq_along(weights)] + if(length(weights) > 2 || anyNA(names(weights)) || !all(names(weights) %in% c("accuracy", "cost")) || anyDuplicated(names(weights))) + stop("'weights' must have one element named \"accuracy\" and one named \"cost\".") + if(abs(sum(weights) - 1) > sqrt(.Machine$double.eps)) + stop("'weights' must sum to 1.") + weightOf <- function(name) if(name %in% names(weights)) weights[[name]] else 0 + + hasCost <- "cost" %in% names(object[["performance"]]) summaryTable <- data.frame(Pathway = rownames(object[["performance"]]), `Balanced Accuracy` = object[["performance"]][, "accuracy"], check.names = FALSE) if(hasCost) summaryTable[, "Total Cost"] <- object[["performance"]][, "cost"] - + + finalScores <- rank(object[["performance"]][, "accuracy"], na.last = "keep") * weightOf("accuracy") if(hasCost) - rankingScores <- list(rank(object[["performance"]][, "accuracy"]), rank(-object[["performance"]][, "cost"])) - else - rankingScores <- list(rank(object[["performance"]][, "accuracy"])) - finalScores <- mapply(function(scores, weight) scores * weight, rankingScores, as.list(weights)) - if(!is(finalScores, "tabular")) finalScores <- matrix(finalScores, nrow = 1) - finalScores <- rowSums(finalScores) + finalScores <- finalScores + rank(-object[["performance"]][, "cost"], na.last = "keep") * weightOf("cost") summaryTable <- cbind(summaryTable, Score = finalScores) summaryTable <- summaryTable[order(summaryTable[, "Score"], decreasing = TRUE), ] summaryTable @@ -396,15 +464,16 @@ bubblePlot <- function (precisionPathways, ...) { #' @export bubblePlot.PrecisionPathways <- function(precisionPathways, pathwayColours = NULL, ...) { - ggplot2::theme_set(ggplot2::theme_classic() + ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA))) if(is.null(pathwayColours)) pathwayColours <- scales::hue_pal()(length(precisionPathways[["pathways"]])) performance <- precisionPathways[["performance"]] performance <- cbind(Sequence = rownames(performance), performance) ggplot2::ggplot(performance, aes(x = accuracy, y = cost, colour = Sequence, size = 4)) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = pathwayColours) + ggplot2::labs(x = "Balanced Accuracy", y = "Total Cost") + - ggplot2::guides(size = "none") + ggplot2::scale_x_continuous(limits = c(0.5, NA)) + ggplot2::scale_y_continuous(limits = c(0, NA)) + ggplot2::guides(size = "none") + ggplot2::scale_y_continuous(limits = c(0, NA)) + + ggplot2::theme_classic() + ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA)) } + #' @export #' @rdname precisionPathwaysEvaluations flowchart <- function (precisionPathways, ...) { @@ -488,6 +557,7 @@ flowchart.PrecisionPathways <- function(precisionPathways, pathway, } } + #' @export #' @rdname precisionPathwaysEvaluations strataPlot <- function (precisionPathways, ...) { @@ -517,18 +587,19 @@ strataPlot.PrecisionPathways <- function(precisionPathways, pathway, classColour samplesTiers <- dplyr::arrange(as.data.frame(samplesTiers), Tier, trueClass, Accuracy) samplesTiers$ID = 1:nrow(samplesTiers) - samplesTiers$colour = ifelse(samplesTiers$trueClass == levels(samplesTiers[, "Predicted"])[1], classColours["class1"], classColours["class2"]) samplesTiers$Tier <- droplevels(samplesTiers$Tier) + # The true classes are drawn as a strip of tiles above the tiers. + samplesTiers$trueClassRow <- length(levels(samplesTiers[, "Tier"])) + 0.8 + trueClassColours <- setNames(unname(classColours[c("class1", "class2")]), levels(samplesTiers[, "Predicted"])[1:2]) strataPlot <- ggplot2::ggplot(mapping = ggplot2::aes(x = ID, y = Tier), data = samplesTiers) + - ggplot2::geom_tile(ggplot2::aes(fill = trueClass)) + - ggplot2::scale_fill_manual(values = unname(classColours)) + - ggplot2::labs(title = paste("Pathway:", pathway), fill = "True Class", x = "", y = "") + - ggplot2::guides(fill = ggplot2::guide_legend(title.position = "top")) + - ggnewscale::new_scale_fill() + ggplot2::geom_tile(ggplot2::aes(fill = Accuracy)) + - ggplot2::scale_fill_gradient2(low = "#86C57C", mid = "white", high = "#FFFF28", midpoint = 0.5) + - ggplot2::labs(fill = "Accuracy") + - ggplot2::guides(fill = ggplot2::guide_colorbar(title.position = "top")) + + ggplot2::scale_fill_gradient2(low = "#86C57C", mid = "white", high = "#FFFF28", midpoint = 0.5, + name = "Accuracy", guide = ggplot2::guide_colorbar(title.position = "top")) + + ggplot2::labs(title = paste("Pathway:", pathway), x = "", y = "") + + ggnewscale::new_scale_fill() + + ggplot2::geom_tile(ggplot2::aes(y = trueClassRow, fill = trueClass), height = 0.6) + + ggplot2::scale_fill_manual(values = trueClassColours, na.value = "grey70", + name = "True Class", guide = ggplot2::guide_legend(title.position = "top")) + ggplot2::theme(panel.background = ggplot2::element_blank(), axis.text.x = ggplot2::element_blank(), aspect.ratio = 1/4, @@ -537,12 +608,7 @@ strataPlot.PrecisionPathways <- function(precisionPathways, pathway, classColour legend.text = ggplot2::element_text(size = 10), legend.position = "bottom", axis.text = ggplot2::element_text(size = 15)) + - ggplot2::annotate("tile", - x = samplesTiers$ID, - y = length(levels(samplesTiers[, "Tier"])) + 0.8, - height = 0.6, - fill = samplesTiers$colour) + ggplot2::coord_cartesian(expand = FALSE) strataPlot -} \ No newline at end of file +} diff --git a/R/utilities.R b/R/utilities.R index 6d182b7..284453e 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -170,12 +170,12 @@ splitsTestInfo <- function(samplesSplits = c("k-Fold", "Permute k-Fold", "Permut tuneCombo <- tuneCombosSelect[rowIndex, , drop = FALSE] if(!identical(names(tuneCombo), "None")) # Add real parameters before function call. paramList <- append(paramList, tuneCombo) - if(attr(featureRanking, "name") == "randomSelection") + if(identical(attr(featureRanking, "name"), "randomSelection")) paramList <- append(paramList, list(nFeatures = topNfeatures)) do.call(featureRanking, paramList) }) - if(attr(featureRanking, "name") %in% c("randomSelection", "previousSelection", "Union Selection")) # Actually selection not ranking. + if(isTRUE(attr(featureRanking, "name") %in% c("randomSelection", "previousSelection", "Union Selection"))) # Actually selection not ranking. return(list(NULL, rankings[[1]], NULL)) if(crossValParams@tuneMode == "none") # No parameters to choose between. @@ -728,10 +728,11 @@ predict.dlda <- function(object, newdata, ...) { # Remove once sparsediscrim is if(!is.null(fixed)) { + if(!is.matrix(permutations)) permutations <- matrix(permutations, ncol = 1) for(rowIndex in seq_len(nrow(fixed))) { keepColumns <- permutations[fixed[rowIndex, 1], ] == fixed[rowIndex, 2] - permutations <- permutations[, keepColumns] + permutations <- permutations[, keepColumns, drop = FALSE] } } permutations diff --git a/tests/testthat/test-classes.R b/tests/testthat/test-classes.R new file mode 100644 index 0000000..3a1e4e0 --- /dev/null +++ b/tests/testthat/test-classes.R @@ -0,0 +1,51 @@ +# Parameter classes: constructors, informative errors and show methods. + +test_that("one-step classifiers give an informative error for PredictParams", { + expect_error(PredictParams("naiveBayes"), "trains and predicts in one function") + expect_error(PredictParams("notAClassifier"), "not a classifier keyword") + expect_output(show(PredictParams("DLDA")), "An object of class 'PredictParams'.", fixed = TRUE) +}) + +test_that("SelectParams accepts 'none' and plain functions", { + expect_null(SelectParams("none")) + expect_s4_class(ModellingParams(selectParams = SelectParams("none")), "ModellingParams") + + fromPackage <- SelectParams(ClassifyR:::differentMeansRanking) + expect_s4_class(fromPackage, "SelectParams") + expect_output(show(fromPackage), "Selection Name: Difference in Means.") + + varianceRanking <- function(measurementsTrain, classesTrain, verbose = 3) + colnames(measurementsTrain)[order(apply(measurementsTrain, 2, var), decreasing = TRUE)] + userParams <- SelectParams(varianceRanking, nFeatures = 5) + expect_output(show(userParams), "User-specified Ranking") + data <- makeTwoClass() + result <- runTests(data$measurements, data$classes, CrossValParams(permutations = 1, folds = 3, parallelParams = BiocParallel::SerialParam()), + ModellingParams(balancing = "none", selectParams = userParams), verbose = 0) + expect_true(all(lengths(chosenFeatureNames(result)) == 5)) +}) + +test_that("balancing is not accepted by TrainParams", { + expect_error(TrainParams("DLDA", balancing = "none"), "ModellingParams") + expect_s4_class(TrainParams("DLDA"), "TrainParams") +}) + +test_that("show methods print complete lines", { + ensemble <- SelectParams(list("t-test", "limma")) + expect_output(show(ensemble), "Minimum Functions Selected By: 1.", fixed = TRUE) + noted <- SelectParams("t-test", characteristics = S4Vectors::DataFrame(characteristic = "Note", value = "x")) + expect_output(show(noted), "Selection Name: Difference in Means.\nNote: x.", fixed = TRUE) + expect_output(show(TrainParams("DLDA")), "Classifier Name: Diagonal LDA.") + expect_output(show(TransformParams("diffLoc")), "Transform Name: Subtraction From Training Set Location.") + sets <- FeatureSetCollection(setNames(lapply(1:5, function(index) paste0("g", 1:8)), paste0("set", 1:5))) + outputFile <- tempfile() + sink(outputFile) + show(sets) + sink() + printed <- readChar(outputFile, file.size(outputFile)) + expect_true(endsWith(printed, "set5: g1, g2, g3, g4, g5, ...\n")) +}) + +test_that("performance metrics record whether higher is better", { + table <- ClassifyR:::.ClassifyRenvir[["performanceInfoTable"]] + expect_true(all(table[, "better"] %in% c("higher", "lower"))) +}) diff --git a/tests/testthat/test-crisscross.R b/tests/testthat/test-crisscross.R new file mode 100644 index 0000000..216c048 --- /dev/null +++ b/tests/testthat/test-crisscross.R @@ -0,0 +1,70 @@ +# Cross-cohort evaluation with crissCrossValidate and its plot. + +makeCohort <- function(seed, nSamples = 40, features = paste0("g", 1:60), shift = 0) +{ + set.seed(seed) + classes <- factor(rep(c("A", "B"), length.out = nSamples)) + measurements <- matrix(rnorm(nSamples * length(features)), nSamples, + dimnames = list(paste0("c", seed, "s", seq_len(nSamples)), features)) + measurements[classes == "B", 1:3] <- measurements[classes == "B", 1:3] + shift + list(measurements = measurements, classes = classes) +} + +test_that("only features measured in every data set are used", { + first <- makeCohort(1, features = paste0("g", 1:60), shift = 1.5) + second <- makeCohort(2, features = paste0("g", c(1:50, 61:70)), shift = 1.5) + set.seed(1) + expect_message(result <- crissCrossValidate(list(first = first$measurements, second = second$measurements), + list(first = first$classes, second = second$classes), + nFeatures = 5, classifier = "DLDA", trainType = "modelTrain"), + "50 features present in every data set") + expect_equal(dim(result[["real"]]), c(2, 2)) + expect_equal(result[["params"]][["diagonal"]], "resubstitution") +}) + +test_that("the modelTest diagonal is cross-validation within the data set", { + # Pure noise: features chosen with the test samples would look predictive. + first <- makeCohort(1, nSamples = 60, features = paste0("g", 1:500)) + second <- makeCohort(2, nSamples = 60, features = paste0("g", 1:500)) + set.seed(1) + result <- crissCrossValidate(list(first = first$measurements, second = second$measurements), + list(first = first$classes, second = second$classes), + nFeatures = 20, classifier = "DLDA", trainType = "modelTest", nRepeats = 3) + expect_equal(result[["params"]][["diagonal"]], "cross-validation") + expect_true(all(diag(result[["real"]]) < 0.7)) + + set.seed(1) + crossValidated <- crossValidate(first$measurements, first$classes, nFeatures = 20, selectionMethod = "t-test", + classifier = "DLDA", multiViewMethod = "none", nFolds = 5, nRepeats = 3) + expect_equal(result[["real"]][1, 1], + round(mean(performance(calcCVperformance(crossValidated, "Balanced Accuracy"))[["Balanced Accuracy"]]), 2)) +}) + +test_that("performance type and outcome type are resolved for every data set", { + first <- makeCohort(1, shift = 1.5) + survival <- survival::Surv(rexp(40), rep(1, 40)) + expect_error(crissCrossValidate(list(first = first$measurements, second = first$measurements), + list(first = first$classes, second = survival), classifier = "DLDA"), + "same kind") +}) + +test_that("random features work with several values of nFeatures, and the plot is returned", { + first <- makeCohort(1, shift = 1.5) + second <- makeCohort(2, shift = 1.5) + set.seed(1) + result <- suppressWarnings(crissCrossValidate(list(first = first$measurements, second = second$measurements), + list(first = first$classes, second = second$classes), + nFeatures = c(5, 10), classifier = "DLDA", trainType = "modelTest", + nRepeats = 2, doRandomFeatures = TRUE)) + expect_equal(dim(result[["random"]]), c(2, 2)) + + searchBefore <- search() + plotted <- crissCrossPlot(result, includeValues = TRUE) + expect_identical(search(), searchBefore) + expect_s3_class(plotted, "ggplot") + + result[["random"]] <- NULL + plotted <- crissCrossPlot(result) + expect_s3_class(plotted, "ggplot") + expect_match(plotted[["labels"]][["caption"]], "cross-validation") +}) diff --git a/tests/testthat/test-pathways.R b/tests/testthat/test-pathways.R new file mode 100644 index 0000000..d3cf9c0 --- /dev/null +++ b/tests/testthat/test-pathways.R @@ -0,0 +1,169 @@ +# Precision pathways: training, prediction of new samples, costs, summaries and plots. + +makePathwaysData <- function(seed, nSamples, assays = c("RNA", "protein")) +{ + set.seed(seed) + outcome <- factor(rep(c("Good", "Poor"), each = nSamples / 2)) + clinical <- S4Vectors::DataFrame(Outcome = outcome, Age = rnorm(nSamples) + 0.8 * (outcome == "Poor"), + Stage = rnorm(nSamples), row.names = paste0("S", seed, "_", seq_len(nSamples))) + experiments <- lapply(seq_along(assays), function(index) + { + values <- matrix(rnorm(30 * nSamples), 30, dimnames = list(paste0(assays[index], "f", 1:30), rownames(clinical))) + values[1:4, outcome == "Poor"] <- values[1:4, outcome == "Poor"] + 1.2 + values + }) + names(experiments) <- assays + MultiAssayExperiment::MultiAssayExperiment(experiments, colData = clinical) +} +useClinical <- list(clinical = c("Age", "Stage")) + +test_that(".permutations keeps a matrix when one permutation remains", { + permutations <- ClassifyR:::.permutations(c("clinical", "RNA"), fixed = data.frame(1, "clinical")) + expect_true(is.matrix(permutations)) + expect_equal(permutations[, 1], c("clinical", "RNA")) + permutations <- ClassifyR:::.permutations(c("clinical", "RNA", "protein"), fixed = data.frame(1, "clinical")) + expect_equal(dim(permutations), c(3, 2)) +}) + +test_that("the default mode works for clinical data and one assay", { + data <- makePathwaysData(1, 40, "RNA") + set.seed(1) + pathways <- precisionPathwaysTrain(data, "Outcome", useFeatures = useClinical, nRepeats = 2, nFolds = 3, + minAssaySamples = 2, classifier = c(clinical = "DLDA", RNA = "DLDA"), + selectionMethod = c(clinical = "none", RNA = "t-test")) + expect_s3_class(pathways, "PrecisionPathways") + expect_true(all(startsWith(names(pathways[["pathways"]]), "clinical"))) +}) + +test_that("fixedAssays may be NULL or several assays, and lists without useFeatures are accepted", { + data <- makePathwaysData(1, 40) + set.seed(1) + noFixed <- precisionPathwaysTrain(data, "Outcome", fixedAssays = NULL, useFeatures = useClinical, nRepeats = 2, + nFolds = 3, classifier = c(clinical = "DLDA", RNA = "DLDA", protein = "DLDA")) + expect_equal(length(noFixed[["assaysPermutations"]]), 6) + set.seed(1) + twoFixed <- precisionPathwaysTrain(data, "Outcome", fixedAssays = c("clinical", "protein"), useFeatures = useClinical, + nRepeats = 2, nFolds = 3, classifier = c(clinical = "DLDA", RNA = "DLDA", protein = "DLDA")) + expect_equal(twoFixed[["assaysPermutations"]], list(c("clinical", "protein", "RNA"))) + + clinical <- as.data.frame(MultiAssayExperiment::colData(data)) + dataList <- list(clinical = clinical[, c("Age", "Stage")], RNA = t(MultiAssayExperiment::assay(data, "RNA"))) + set.seed(1) + expect_warning(fromList <- precisionPathwaysTrain(dataList, clinical[, "Outcome"], nRepeats = 2, nFolds = 3, + classifier = c(clinical = "DLDA", RNA = "DLDA")), "useFeatures") + expect_s3_class(fromList, "PrecisionPathways") + expect_equal(fromList[["useFeatures"]][["clinical"]], c("Age", "Stage")) +}) + +trainSmall <- function(classifier, data = makePathwaysData(1, 60)) +{ + set.seed(2) + precisionPathwaysTrain(data, "Outcome", useFeatures = useClinical, nRepeats = 3, nFolds = 3, minAssaySamples = 3, + mode = "combinatorial", classifier = classifier, + selectionMethod = c(clinical = "none", RNA = "t-test", protein = "t-test")) +} + +test_that("models are matched to classifiers and pathways by assay name, not by position", { + newData <- makePathwaysData(2, 30) + inOrder <- trainSmall(c(clinical = "GLM", protein = "randomForest", RNA = "DLDA")) + expect_equal(names(inOrder[["models"]]), c("clinical", "RNA", "protein")) + reordered <- inOrder + reordered[["parameters"]][["classifier"]] <- inOrder[["parameters"]][["classifier"]][c("RNA", "protein", "clinical")] + reordered[["models"]] <- inOrder[["models"]][c("protein", "clinical", "RNA")] + set.seed(5) # Random forest breaks tied votes at random. + predictedInOrder <- precisionPathwaysPredict(inOrder, newData, "Outcome") + set.seed(5) + predictedReordered <- precisionPathwaysPredict(reordered, newData, "Outcome") + expect_equal(lapply(predictedInOrder[["pathways"]], "[[", "individuals"), + lapply(predictedReordered[["pathways"]], "[[", "individuals")) + + # Each new sample is predicted by every fold model of every repeat. + individuals <- predictedInOrder[["pathways"]][[1]][["individuals"]] + expect_setequal(individuals[, "Sample ID"], rownames(MultiAssayExperiment::colData(newData))) + # Fold models keep only what prediction needs. + expect_null(attr(models(inOrder[["models"]][["protein"]])[[1]], "forImportance")) +}) + +test_that("new samples can be predicted without their classes", { + pathways <- trainSmall(c(clinical = "DLDA", RNA = "DLDA", protein = "DLDA")) + newData <- makePathwaysData(2, 30) + withClasses <- precisionPathwaysPredict(pathways, newData, "Outcome") + withoutClasses <- precisionPathwaysPredict(pathways, newData) + expect_equal(lapply(withoutClasses[["pathways"]], function(pathway) pathway[["individuals"]][, "Predicted"]), + lapply(withClasses[["pathways"]], function(pathway) pathway[["individuals"]][, "Predicted"])) + costed <- calcCostsAndPerformance(withoutClasses, c(clinical = 50, RNA = 700, protein = 400)) + expect_true(all(is.na(costed[["performance"]][, "accuracy"]))) + expect_true(all(costed[["performance"]][, "cost"] > 0)) +}) + +test_that("assay names containing '-' are kept intact", { + data <- makePathwaysData(1, 40, c("RNA-seq", "protein")) + set.seed(1) + pathways <- suppressWarnings(precisionPathwaysTrain(data, "Outcome", useFeatures = useClinical, nRepeats = 2, nFolds = 3, + minAssaySamples = 3, classifier = setNames(rep("DLDA", 3), c("clinical", "RNA-seq", "protein")), + selectionMethod = c(clinical = "none", `RNA-seq` = "t-test", protein = "t-test"))) + expect_true(all(unlist(lapply(pathways[["pathways"]], "[[", "assays")) %in% c("clinical", "RNA-seq", "protein"))) + predicted <- suppressWarnings(precisionPathwaysPredict(pathways, makePathwaysData(2, 30, c("RNA-seq", "protein")), "Outcome")) + expect_s3_class(predicted, "PrecisionPathways") +}) + +test_that("the confidence cut-off is a minimum", { + counts <- as.table(matrix(c(18, 2, 2, 18, 17, 3), ncol = 2, byrow = TRUE, + dimnames = list(c("a", "b", "c"), c("Good", "Poor")))) + expect_equal(ClassifyR:::.pathwaysConfident(counts, 0.8), c("a", "b")) +}) + +# A pathway built by hand: two samples predicted by clinical data and one passed on to RNA. +handPathways <- function(withTested = TRUE) +{ + individuals <- S4Vectors::DataFrame(Tier = factor(c("clinical", "clinical", "RNA"), levels = c("clinical", "RNA")), + `Sample ID` = c("s1", "s2", "s3"), + Predicted = factor(c("A", "B", "B"), levels = c("A", "B")), + Accuracy = c(1, 1, 1), check.names = FALSE) + pathway <- list(pathway = "clinical-RNA", assays = c("clinical", "RNA"), assaysTested = c("clinical", "RNA"), + individuals = individuals, + tiers = S4Vectors::DataFrame(Tier = c("clinical", "RNA"), `Balanced Accuracy` = c(1, 1), check.names = FALSE)) + if(!withTested) pathway[c("assays", "assaysTested")] <- NULL + structure(list(pathways = list(`clinical-RNA` = pathway), + testSampleInfo = S4Vectors::DataFrame(sampleID = c("s1", "s2", "s3"), class = factor(c("A", "B", "B")))), + class = "PrecisionPathways") +} + +test_that("a sample is charged for every tier it passed through", { + costs <- c(clinical = 50, RNA = 700) + expect_equal(calcCostsAndPerformance(handPathways(), costs)[["performance"]][, "cost"], 50 + 50 + (50 + 700)) + expect_equal(calcCostsAndPerformance(handPathways(FALSE), costs)[["performance"]][, "cost"], 850) + expect_equal(calcCostsAndPerformance(handPathways(), c(RNA = 700))[["performance"]][, "cost"], 700) +}) + +test_that("summary weights are matched by name and must sum to 1", { + pathways <- structure(list(performance = data.frame(accuracy = c(0.9, 0.6, 0.55), cost = c(9000, 100, 50), + row.names = c("expensiveGood", "cheapA", "cheapest"))), + class = "PrecisionPathways") + costFirst <- summary(pathways, weights = c(cost = 0.9, accuracy = 0.1)) + accuracyFirst <- summary(pathways, weights = c(accuracy = 0.1, cost = 0.9)) + expect_equal(costFirst, accuracyFirst) + expect_equal(costFirst[1, "Pathway"], "cheapest") + expect_equal(summary(pathways, weights = c(accuracy = 1, cost = 0))[1, "Pathway"], "expensiveGood") + expect_error(summary(pathways, weights = c(accuracy = 5, cost = 5)), "sum to 1") + expect_error(summary(pathways, weights = c(accuracy = 0.5, price = 0.5)), "named") +}) + +test_that("plots keep the global theme, every pathway and the true classes", { + pathways <- handPathways() + pathways[["performance"]] <- data.frame(accuracy = 0.3, cost = 850, row.names = "clinical-RNA") + pathways[["models"]] <- list(ClassifyResult(S4Vectors::DataFrame(characteristic = "Classifier Name", value = "DLDA"), + c("s1", "s2", "s3"), "f1", list(), list(), list(), NULL, + S4Vectors::DataFrame(sample = c("s1", "s2", "s3"), class = factor(c("A", "B", "B"))), + factor(c("A", "B", "B")))) + themeBefore <- ggplot2::theme_get() + bubbles <- bubblePlot(pathways) + expect_identical(ggplot2::theme_get(), themeBefore) + expect_equal(nrow(ggplot2::layer_data(bubbles, 1)), 1) # Accuracy below 0.5 is still drawn. + + strata <- strataPlot(pathways, "clinical-RNA") + expect_identical(ggplot2::theme_get(), themeBefore) + trueClassTiles <- ggplot2::layer_data(strata, 2) + expect_equal(nrow(trueClassTiles), 3) + expect_true(all(trueClassTiles[, "y"] > 2)) # Drawn above the tiers, not underneath them. +})