10/ Metrics and plots: no internal rounding, tied risks count 1/2, fold grouping for all metrics, chance line per metric (changes results slightly)

R/ROCplot.R

@@ -203,10 +203,10 @@ setMethod("ROCplot", "list",
     if(numberDistinctClasses == 2)
     {
       plotData <- subset(plotData, class == distinctClasses[2])
-      plotData[, comparisonName] <- paste(plotData[, comparisonName], " (AUC ", plotData[, "AUC"], ')', sep = '')
+      plotData[, comparisonName] <- paste(plotData[, comparisonName], " (AUC ", round(plotData[, "AUC"], 2), ')', sep = '')
       plotData[, comparisonName] <- factor(plotData[, comparisonName], levels = unique(plotData[, comparisonName]))
     } else {
-      plotData[, "class"] <- paste(plotData[, "class"], " (AUC ", plotData[, "AUC"], ')', sep = '')
+      plotData[, "class"] <- paste(plotData[, "class"], " (AUC ", round(plotData[, "AUC"], 2), ')', sep = '')
       plotData[, "class"] <- factor(plotData[, "class"], levels = unique(plotData[, "class"]))
     }
   }
@@ -231,7 +231,7 @@ setMethod("ROCplot", "list",
                   ROCplot <- ROCplot + ggplot2::xlab(xLabel) + ggplot2::ylab(yLabel) + ggplot2::ggtitle(plotTitle) + ggplot2::theme(axis.title = ggplot2::element_text(size = fontSizes[2]), plot.title = ggplot2::element_text(size = fontSizes[1], hjust = 0.5))
 
                 if(length(results) == 1 && showAUC == TRUE && numberDistinctClasses == 2)
-                  ROCplot <- ROCplot + ggplot2::annotate("text", x = Inf, y = 0, label = paste("AUC =", plotData[1, "AUC"]), hjust = 1.1, size = fontSizes[2] * 5/14) + ggplot2::theme(legend.position = "none")
+                  ROCplot <- ROCplot + ggplot2::annotate("text", x = Inf, y = 0, label = paste("AUC =", round(plotData[1, "AUC"], 2)), hjust = 1.1, size = fontSizes[2] * 5/14) + ggplot2::theme(legend.position = "none")
 
                 if(length(results) > 1 && numberDistinctClasses > 2)
                   ROCplot <- ROCplot + ggplot2::facet_wrap(ggplot2::vars(!!comparison))

R/calcPerformance.R

@@ -50,25 +50,30 @@
 #' \item{\code{"Micro F1"}: F1 score obtained by calculating the
 #' harmonic mean of micro precision and micro recall.}
 #' \item{\code{"Macro Precision"}: Sum of the ratios of the number of correct predictions in each
-#' class to the number of samples predicted to be in each class, divided by the number of classes.}
+#' class to the number of samples predicted to be in each class, divided by the number of classes. A class which is never
+#' predicted has a precision of 0.}
 #' \item{\code{"Macro Recall"}: Sum of the ratios of the number of correct predictions
 #' in each class to the number of samples in each class, divided by the number of classes.}
 #' \item{\code{"Macro F1"}: F1 score obtained by calculating the harmonic mean of macro precision
 #' and macro recall.}
 #' \item{\code{"Matthews Correlation Coefficient"}: Matthews Correlation Coefficient (MCC). A score
 #' between -1 and 1 indicating how concordant the predicted classes are to the actual classes. Only defined if
 #' there are two classes.}
-#' \item{\code{"AUC"}: Area Under the Curve. An area ranging from 0 to 1, under the ROC.}
+#' \item{\code{"AUC"}: Area Under the Curve. An area ranging from 0 to 1, under the ROC. For more than two classes,
+#'         the mean of the one-versus-rest areas. If a group of predictions lacks a sample of some class, its AUC is
+#'         \code{NaN}, with a warning.}
 #' \item{\code{"C-index"}: For survival data, the concordance index, for models which produce risk scores. Ranges from 0 to 1.}
-#' \item{\code{"Sample C-index"}: Per-individual C-index.}
+#' \item{\code{"Sample C-index"}: Per-individual C-index. For each sample, the proportion of comparable pairs it
+#'         belongs to which are concordant, with a pair of tied risk scores counting as half concordant, as for the C-index.}
 #' }
 #' 
 #' @param actualOutcome A factor vector or survival information specifying each sample's known outcome.
 #' @param predictedOutcome A factor vector or survival information of the same length as \code{actualOutcome} specifying each sample's predicted outcome.
 #' @param grouping Default: \code{"permutation"}. If the cross-validation was k-fold, then this determines whether the metric will be calculated for samples
 #' grouped by permutation or by fold, if the value is \code{"fold"}. For small sample sizes, \code{"permutation"} would suit. But, for large sample sizes,
-#' \code{"fold"} would be preferable, as class membership probabilities or risk scores are not directly comparable between folds. This setting makes
-#' no difference to error or accuracy metrics, apart from their variability.
+#' \code{"fold"} would be preferable, as class membership probabilities or risk scores are not directly comparable between folds. If
+#' \code{"fold"}, the metric is calculated for each fold and the folds' values are averaged, giving one value for each permutation.
+#' Metrics are not rounded.
 #' 
 #' @return If \code{calcCVperformance} was run, an updated
 #' \code{\linkS4class{ClassifyResult}} object, with new metric values in the
@@ -225,6 +230,8 @@ setMethod("calcCVperformance", "ClassifyResult",
         predictedOutcome <- factor(result@predictions[, "class"], levels = classLevels)
         actualOutcome <- factor(actualOutcomeOrdered, levels = classLevels, ordered = TRUE)
         performance <- .calcPerformance(actualOutcome, predictedOutcome, samples, performanceType, groupID)
+        if(!performanceType %in% c("Sample Error", "Sample Accuracy") && grepl(':', names(performance[["values"]])[1])) # Then average for each permutation.
+          performance[["values"]] <- .averagePermutations(performance[["values"]])
         result@performance[[performance[["name"]]]] <- performance[["values"]]
       }
   }
@@ -282,17 +289,22 @@ setMethod("calcCVperformance", "ClassifyResult",
       earlierEvent <- outer(times, times, '<') & iterationSurv[, "status"] == 1
       higherRisk <- outer(iterationPredictions, iterationPredictions, '>')
       lowerRisk <- outer(iterationPredictions, iterationPredictions, '<')
+      tiedRisk <- outer(iterationPredictions, iterationPredictions, '==')
       # Concordant: the sample with the earlier event has the higher risk. Discordant: it has the lower risk.
+      # Tied: both have the same risk, which counts as half concordant, as for the C-index.
       concordant <- earlierEvent & higherRisk
       discordant <- earlierEvent & lowerRisk
+      tied <- earlierEvent & tiedRisk
       data.frame(sample = iterationSamples,
                  concordant = rowSums(concordant) + colSums(concordant),
-                 discordant = rowSums(discordant) + colSums(discordant))
+                 discordant = rowSums(discordant) + colSums(discordant),
+                 tied = rowSums(tied) + colSums(tied))
     }, actualOutcome, predictedOutcome, samples, SIMPLIFY = FALSE))
 
     concordants <- rowsum(pairCounts[, "concordant"], pairCounts[, "sample"])[, 1]
     discordants <- rowsum(pairCounts[, "discordant"], pairCounts[, "sample"])[, 1]
-    Cindex <- round(concordants / (concordants + discordants), 2)[allSamples]
+    tieds <- rowsum(pairCounts[, "tied"], pairCounts[, "sample"])[, 1]
+    Cindex <- ((concordants + 0.5 * tieds) / (concordants + discordants + tieds))[allSamples]
     names(Cindex) <- allSamples
     Cindex[is.nan(Cindex)] <- NA # A censored individual with a time shorter than all events has no comparable pairs.
     return(list(name = performanceType, values = Cindex))
@@ -329,8 +341,13 @@ setMethod("calcCVperformance", "ClassifyResult",
         mean(classErrors / classSizes, na.rm = TRUE)
     }, actualOutcome, predictedOutcome, SIMPLIFY = FALSE))
   } else if(performanceType == "AUC") {
+    missingClass <- sapply(actualOutcome, function(iterationClasses) any(table(iterationClasses) == 0))
+    if(any(missingClass))
+      warning(sum(missingClass), " of ", length(missingClass), " groups of predictions lack a sample of some class. ",
+              "Their AUC is NaN.")
     performanceValues <- unlist(mapply(function(iterationClasses, iterationPredictions)
     {
+      if(any(table(iterationClasses) == 0)) return(NaN)
       classesTable <- do.call(rbind, lapply(levels(iterationClasses), function(class)
       {
         rates <- .ROCrates(iterationPredictions[, class], iterationClasses == class)
@@ -374,7 +391,9 @@ setMethod("calcCVperformance", "ClassifyResult",
       }
       if(performanceType %in% c("Macro Precision", "Macro F1"))
       {
-        macroP <- sum(truePositives / (truePositives + falsePositives)) / nrow(confusionMatrix)
+        classPrecisions <- truePositives / (truePositives + falsePositives)
+        classPrecisions[truePositives + falsePositives == 0] <- 0 # The class was never predicted.
+        macroP <- sum(classPrecisions) / nrow(confusionMatrix)
         if(performanceType == "Macro Precision") return(macroP)
       }
       if(performanceType %in% c("Macro Recall", "Macro F1"))

R/crissCrossValidate.R

@@ -156,7 +156,6 @@ crissCrossValidate <- function(measurements, outcomes,
             dimnames = list(paste("Select and Train", names(measurements)),
                             paste("Predict",         names(measurements)))
         )
-        realPerformance <- round(realPerformance, 2)
 
     } else {
         # trainType == "modelTest"
@@ -208,7 +207,6 @@ crissCrossValidate <- function(measurements, outcomes,
             dimnames = list(paste("Select", names(measurements)),
                             paste("Cross-validate", names(measurements)))
         )
-        realPerformance <- round(realPerformance, 2)
     }
 
     result <- list(real = realPerformance) # Store the results
@@ -244,7 +242,6 @@ crissCrossValidate <- function(measurements, outcomes,
             dimnames = list(paste("Random Select", names(measurements)),
                             paste("Cross-validate", names(measurements)))
         )
-        randomPerformance <- round(randomPerformance, 2)
 
         result$random <- randomPerformance
     }
@@ -294,7 +291,7 @@ crissCrossValidate <- function(measurements, outcomes,
 
         topMatrix <- matrix(topPerformance, nrow = 1,
                             dimnames = list("TOP", names(top_measurements)))
-        result$top <- round(topMatrix, 2)
+        result$top <- topMatrix
     }
 
     # Store the parameters
@@ -330,15 +327,16 @@ crissCrossValidate <- function(measurements, outcomes,
 #' A function to plot the output of the crissCrossValidate function.
 #'
 #' This function generates a heatmap of the cross-validation results from
-#' \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.
+#' \code{\link{crissCrossValidate}}. The diagonal, where the features are chosen or the model is trained in the data set
+#' that is tested, is hidden unless \code{showDiagonal = TRUE}. For \code{trainType = "modelTrain"}, it is
+#' resubstitution. For \code{trainType = "modelTest"}, it is cross-validation within each data set. If shown, it is
+#' outlined in black and a caption says how it was evaluated.
 #'
 #' @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 of a \code{"modelTrain"} result
-#'        are set to \code{NA} and appear grayed-out or blank. Defaults to \code{FALSE}.
+#' @param showDiagonal Logical. If \code{FALSE}, the diagonal cells are set to \code{NA} and appear grayed-out.
+#'        Defaults to \code{FALSE}.
+#' @param showResubMetric Deprecated name of \code{showDiagonal}.
 #' @return A \code{ggplot} object, or a combined plot of two heatmaps if random features were evaluated.
 #'
 #' @import ggplot2
@@ -347,8 +345,14 @@ crissCrossValidate <- function(measurements, outcomes,
 #' @export
 crissCrossPlot <- function(crissCrossResult,
                            includeValues    = FALSE,
-                           showResubMetric  = FALSE)
+                           showDiagonal     = FALSE,
+                           showResubMetric)
 {
+    if(!missing(showResubMetric))
+    {
+        warning("'showResubMetric' is deprecated. Use 'showDiagonal' instead.")
+        showDiagonal <- showResubMetric
+    }
     params <- crissCrossResult[["params"]]
     scalebar_title <- params$performanceType
     # Results created before the diagonal was recorded have a resubstitution diagonal.
@@ -360,7 +364,7 @@ crissCrossPlot <- function(crissCrossResult,
         # Convert to matrix if needed
         mat <- as.matrix(mat)
 
-        if (isResubstitution && !showResubMetric) {
+        if (!showDiagonal) {
             diag(mat) <- NA
         }
 
@@ -389,7 +393,7 @@ crissCrossPlot <- function(crissCrossResult,
             ) +
             coord_fixed()
 
-        if (!isResubstitution || showResubMetric) {
+        if (showDiagonal) {
             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))

R/performancePlot.R

@@ -9,6 +9,11 @@
 #' \code{characteristicsList[["x"]]} is one of \code{aggregate}, its values
 #' are averaged to a single number and a bar is plotted.
 #' 
+#' A dashed line shows the metric's expected value for predictions made by chance: 0.5 for AUC and C-index, 0 for
+#' Matthews Correlation Coefficient and, with k classes, 1/k for accuracy, precision, recall and F1 metrics and 1 - 1/k for
+#' error metrics. This is the value expected from predicting classes uniformly at random. No line is drawn for
+#' other metrics.
+#' 
 #' @aliases performancePlot performancePlot,list-method
 #' @param results A list of \code{\link{ClassifyResult}} objects.
 #' @param aggregate A character vector of the levels of
@@ -145,7 +150,7 @@ setMethod("performancePlot", "list",
     results <- lapply(results, function(result) calcCVperformance(result, metric))
   }
 
-  ifelse(metric == "Matthews Correlation Coefficient", baseline <- 0, baseline <- 0.5)
+  baseline <- .chanceLevel(metric, results[[1]])
 
   plotData <- do.call(rbind, mapply(function(result, index)
                     {
@@ -182,8 +187,8 @@ setMethod("performancePlot", "list",
   legendPosition <- ifelse(showLegend == TRUE, "right", "none")
   characteristicsList <- lapply(characteristicsList, rlang::sym)
 
-  performancePlot <- ggplot2::ggplot() + ggplot2::theme_classic() + ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA)) +
-                     ggplot2::geom_hline(yintercept = baseline, linetype = 2)
+  performancePlot <- ggplot2::ggplot() + ggplot2::theme_classic() + ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA))
+  if(!is.na(baseline)) performancePlot <- performancePlot + ggplot2::geom_hline(yintercept = baseline, linetype = 2)
 
   if(!is.null(yLimits) && rotate90 == FALSE) performancePlot <- performancePlot + ggplot2::coord_cartesian(ylim = yLimits)
   if("fillColour" %in% names(characteristicsList))
@@ -215,3 +220,17 @@ setMethod("performancePlot", "list",
 
   performancePlot
 })
+
+# Expected value of a metric for predictions made by chance. For classes, predictions are uniformly random.
+.chanceLevel <- function(metric, result)
+{
+  if(metric %in% c("AUC", "C-index", "Sample C-index")) return(0.5)
+  if(metric == "Matthews Correlation Coefficient") return(0)
+  outcome <- actualOutcome(result)
+  if(!is.factor(outcome)) return(NA)
+  nClasses <- length(levels(outcome))
+  if(metric %in% c("Accuracy", "Balanced Accuracy", "Sample Accuracy", "Micro Precision", "Micro Recall", "Micro F1",
+                   "Macro Precision", "Macro Recall", "Macro F1")) return(1 / nClasses)
+  if(metric %in% c("Error", "Balanced Error", "Sample Error")) return(1 - 1 / nClasses)
+  NA
+}

R/rankingPlot.R

@@ -9,7 +9,9 @@
 #' overlap between all possible pairs of results. The second kind of summary is
 #' the pair-wise overlap of each level of the comparison factor that is not the
 #' reference level against the reference level. The overlaps are converted to
-#' percentages and plotted as lineplots.
+#' percentages and plotted as lineplots. The overlap of two rankings at a number of top features is the number of features
+#' in both lists of top features, divided by the length of the shorter list (which is shorter than the number of top
+#' features only if a ranking has fewer features than it).
 #' 
 #' If \code{comparison} is \code{"within"}, then the feature selection overlaps
 #' are compared within a particular analysis. The result will inform how stable
@@ -60,7 +62,8 @@
 #' @param xLabelPositions Locations where to put labels on the x-axis.
 #' @param yMax The maximum value of the percentage to plot.
 #' @param title An overall title for the plot.
-#' @param yLabel Label to be used for the y-axis of overlap percentages.
+#' @param yLabel Label to be used for the y-axis of overlap percentages. By default, it states that the percentage of
+#' top features in common is shown.
 #' @param margin The margin to have around the plot.
 #' @param showLegend If \code{TRUE}, a legend is plotted next to the plot. If
 #' FALSE, it is hidden.
@@ -121,7 +124,7 @@ setMethod("rankingPlot", "list",
                                     fonts = c(24, 16, 12, 12, 12, 16)),
                    lineColours = NULL, xLabelPositions = seq(10, 100, 10), yMax = 100,
                    title = if(comparison[1] == "within") "Feature Ranking Stability" else "Feature Ranking Commonality",
-                   yLabel = if(is.null(referenceLevel)) "Average Common Features (%)" else paste("Average Common Features with", referenceLevel, "(%)"),
+                   yLabel = if(is.null(referenceLevel)) "Average Top Features in Common (%)" else paste("Average Top Features in Common with", referenceLevel, "(%)"),
                    margin = grid::unit(c(1, 1, 1, 1), "lines"),
                    showLegend = TRUE, parallelParams = bpparam())
 {
@@ -184,7 +187,7 @@ setMethod("rankingPlot", "list",
         {
           sapply(topRanked, function(top)
           {
-            length(intersect(head(features, top), head(other, top))) / top * 100
+            .topOverlap(features, other, top)
           })
         })
       }, rankedList[1:(length(rankedList) - 1)], 1:(length(rankedList) - 1), SIMPLIFY = FALSE)))
@@ -234,7 +237,7 @@ setMethod("rankingPlot", "list",
             {
               sapply(topRanked, function(top)
               {
-                length(intersect(head(rankings, top), head(otherRanked, top))) / top * 100
+                .topOverlap(rankings, otherRanked, top)
               })          
             })
           })))
@@ -291,3 +294,12 @@ setMethod("rankingPlot", "list",
 
   overlapPlot
 })
+
+# Percentage of the top features of two rankings which are in common. The denominator is the length of the shorter
+# list of top features, so that a ranking with fewer features than top is not penalised.
+.topOverlap <- function(ranking, otherRanking, top)
+{
+  ranking <- head(ranking, top)
+  otherRanking <- head(otherRanking, top)
+  length(intersect(ranking, otherRanking)) / min(length(ranking), length(otherRanking)) * 100
+}

R/runTest.R

@@ -286,7 +286,7 @@ input data. Autmomatically reducing to smaller number.")
 
     if(!is.null(ncol(predictedOutcome)))
         predictedOutcome <- predictedOutcome[, na.omit(match(c("class", "risk"), colnames(predictedOutcome)))]
-    performanceChanges <- round(performancesWithoutEach - calcExternalPerformance(outcomeTest, predictedOutcome, performanceType), 2)
+    performanceChanges <- performancesWithoutEach - calcExternalPerformance(outcomeTest, predictedOutcome, performanceType)
 
     if(is.null(S4Vectors::mcols(measurementsTrain)) || !any(c("assay", "feature") %in% colnames(S4Vectors::mcols(measurementsTrain))))
     {