R/metrics.R

In rtemis: Schema-Driven Machine Learning and Visualization

# metrics.R
# ::rtemis::
# 2019- EDG rtemis.org

#' Error functions
#'
#' Convenience functions for calculating loss.
#' These can be passed as arguments to learners that support custom loss functions.
#'
#' @rdname error
#' @param x Vector of True values
#' @param y Vector of predicted values
#' @param na.rm Logical: If TRUE, remove NA values before computation.
#'
#' @author EDG
#' @keywords internal
#' @noRd
mae <- function(x, y, na.rm = TRUE) {
  error <- x - y
  mean(abs(error), na.rm = na.rm)
} # /rtemis::mae

#' @rdname error
#' @keywords internal
#' @noRd
mse <- function(x, y, na.rm = TRUE) {
  error <- x - y
  mean(error^2, na.rm = na.rm)
} # /rtemis::mse

#' Weighted MSE
#'
#' @rdname error
#' @keywords internal
#' @noRd
msew <- function(x, y, weights = rep(1, length(y)), na.rm = TRUE) {
  error <- x - y
  error <- error * weights
  mean(error^2, na.rm = na.rm)
} # /rtemis::msew

#' @rdname error
#' @keywords internal
#' @noRd
rmse <- function(x, y, na.rm = TRUE) {
  sqrt(mse(x, y, na.rm = na.rm))
} # /rtemis::rmse

#' R-squared
#'
#' @param x Float, vector: True values
#' @param y Float, vector: predicted values
#' @author EDG
#' @keywords internal
#' @noRd
rsq <- function(x, y) {
  SSE <- sum((x - y)^2)
  # Sum of Squares due to Regression (SSR) a.k.a. Explained Sum of Squares (ESS)
  # SSR <- sum((mean(x) - y)^2)
  # Total Sum of Squares (TSS or SST)
  SST <- sum((x - mean(x))^2)
  # R-squared a.k.a. Coefficient of Determination i.e. percent variance explained
  1 - (SSE / SST)
} # /rtemis::rsq


#' Log Loss for a binary classifier
#'
#' @param true_int Integer vector, {0, 1}: True labels (1 is the positive class).
#' @param predicted_prob Float, vector: predicted probabilities.
#' @param eps Float: Small value to prevent log(0).
#'
#' @author EDG
#' @keywords internal
#' @noRd
logloss <- function(true_int, predicted_prob, eps = 1e-16) {
  predicted_prob <- pmax(pmin(predicted_prob, 1 - eps), eps)
  -mean(
    true_int * log(predicted_prob) + (1 - true_int) * log(1 - predicted_prob)
  )
} # /rtemis::logloss

#' Factor harmonize
#'
#' @param reference Reference factor.
#' @param x Input factor.
#' @param verbosity Integer: Verbosity level.
#'
#' @return Factor: x with levels in the same order as reference.
#'
#' @author EDG
#'
#' @keywords internal
#' @noRd
factor_harmonize <- function(reference, x, verbosity = 1L) {
  if (!is.factor(x) || !is.factor(reference)) {
    rtemis.core::abort(
      "Inputs must be factors.",
      class = c("rtemis_type_error", "rtemis_input_error")
    )
  }
  if (!all(levels(x) == levels(reference))) {
    if (!all(levels(x) %in% levels(reference))) {
      if (verbosity > 0L) {
        msg("Levels of x:")
      }
      levels(x)
      if (verbosity > 0L) {
        msg("levels of reference:")
      }
      levels(reference)
      rtemis.core::abort(
        "Levels of two inputs do not match.",
        class = c("rtemis_level_mismatch", "rtemis_data_error")
      )
    }
    if (verbosity > 0L) {
      msg("Input factor levels are not in the same order, correcting")
    }
    x <- factor(x, levels = levels(reference))
  }
  x
} # /rtemis::factor_harmonize


#' F1 score
#'
#' Calculate the F1 score for classification:
#'
#' \deqn{F1 = 2 \frac{Recall \cdot Precision}{Recall + Precision}}{F1 = 2 * (Recall * Precision)/(Recall + Precision)}
#'
#' @param recall Float \[0, 1\]: Recall a.k.a. Sensitivity
#' @param precision Float \[0, 1\]: Precision a.k.a. Positive Predictive Value
#'
#' @author EDG
#' @keywords internal
#' @noRd
f1 <- function(precision, recall) {
  2 * (recall * precision) / (recall + precision)
} # /rtemis::f1


# auc.R
# ::rtemis::
# 2019-23 EDG rtemis.org

#' Area under the ROC Curve
#'
#' Get the Area under the ROC curve to assess classifier performance.
#'
#' @param true_int Integer vector: True labels of outcomes (e.g. c(0, 1, 1))
#' @param predicted_prob Numeric Vector: Probabilities or model scores
#' (e.g. c(.32, .75, .63), etc)
#' @param method Character {"lightAUC" "ROCR"}: Package to use.
#' @param verbosity Integer: Verbosity level.
#'
#' @return Numeric.
#'
#' @author EDG
#' @keywords internal
#' @noRd
#'
#' @examples
#' preds <- c(0.7, 0.55, 0.45, 0.25, 0.6, 0.7, 0.2)
#' labels <- 2L - as.integer(factor(c("a", "a", "a", "b", "b", "b", "b")))
#' auc(labels, preds, method = "lightAUC")
#' auc(labels, preds, method = "ROCR")
auc <- function(
  true_int,
  predicted_prob,
  method = c("lightAUC", "ROCR"),
  verbosity = 0L
) {
  # Checks ----
  method <- match.arg(method)
  check_inherits(true_int, "integer")
  check_float01inc(predicted_prob)
  # method <- match.arg(method)
  if (length(unique(true_int)) == 1) {
    return(NaN)
  }
  # Resolve backend ----
  # AUC is one optional metric computed via an optional (Suggests) backend.
  # If the requested backend cannot be loaded (e.g. lightAUC's compiled stack
  # fails to load after an upstream RcppParallel ABI bump), fall back rather
  # than abort: an unavailable metric backend must not crash model training.
  if (method == "lightAUC" && !auc_backend_available("lightAUC")) {
    if (verbosity > 0L) {
      msg("lightAUC unavailable; falling back to ROCR for AUC.")
    }
    method <- "ROCR"
  }
  if (method == "ROCR" && !auc_backend_available("ROCR")) {
    if (verbosity > 0L) {
      msg("No AUC backend available (lightAUC, ROCR); returning NA for AUC.")
    }
    return(NaN)
  }

  auc. <- if (method == "lightAUC") {
    try(lightAUC::lightAUC(probs = predicted_prob, actuals = true_int))
  } else {
    .pred <- try(ROCR::prediction(
      predicted_prob,
      true_int,
      label.ordering = NULL
    ))
    try(ROCR::performance(.pred, "auc")@y.values[[1]])
  }

  if (inherits(auc., "try-error")) {
    auc. <- NaN
  }

  if (verbosity > 0L) {
    msg("AUC =", auc.)
  }
  auc.
} # /rtemis::auc


# %% auc_backend_available ----
#' Check whether an AUC backend package can be loaded
#'
#' Thin wrapper around `requireNamespace()`. Kept as a standalone internal so
#' `auc()`'s fallback path can be exercised in tests via mocked bindings.
#'
#' @param pkg Character: Package name.
#'
#' @return Logical: `TRUE` if the namespace can be loaded.
#'
#' @author EDG
#' @keywords internal
#' @noRd
auc_backend_available <- function(pkg) {
  requireNamespace(pkg, quietly = TRUE)
} # /rtemis::auc_backend_available


#' Area under the Curve by pairwise concordance
#'
#' Get the Area under the ROC curve to assess classifier performance using pairwise concordance
#'
#' The first level of `true.labels` must be the positive class, and high numbers in
#' `estimated.score` should correspond to the positive class.
#'
#' @param estimated.score Float, Vector: Probabilities or model scores
#' (e.g. c(.32, .75, .63), etc)
#' @param true.labels True labels of outcomes (e.g. c(0, 1, 1))
#' @param verbosity Integer: Verbosity level.
#'
#' @examples
#' true.labels <- factor(c("a", "a", "a", "b", "b", "b", "b"))
#' estimated.score <- c(0.7, 0.55, 0.45, 0.25, 0.6, 0.7, 0.2)
#' auc_pairs(estimated.score, true.labels, verbosity = 1L)
#'
#' @keywords internal
#' @noRd
auc_pairs <- function(estimated.score, true.labels, verbosity = 1L) {
  true.labels <- as.factor(true.labels)
  true.levels <- levels(true.labels)
  n.levels <- length(true.levels)
  if (n.levels == 2) {
    outer.diff <- outer(
      estimated.score[true.labels == true.levels[1]],
      estimated.score[true.labels == true.levels[2]],
      "-"
    )
    .auc <- mean((outer.diff > 0) + .5 * (outer.diff == 0))
  } else {
    rtemis.core::abort(
      "Multiclass AUC does not have a unique definition and is not yet implemented.",
      class = "rtemis_unsupported_error"
    )
  }
  if (verbosity > 0L) {
    msg("Positive class:", true.levels[1])
    msg("AUC =", .auc)
  }
  invisible(.auc)
} # /rtemis::auc_pairs


#' Brier_Score
#'
#' Calculate the Brier_Score for classification:
#'
#' \deqn{BS = \frac{1}{N} \sum_{i=1}^{N} (y_i - p_i)^2}{BS = 1/N * sum_{i=1}^{N} (y_i - p_i)^2}
#'
#' @param true_int Integer vector, {0, 1}: True labels
#' @param predicted_prob Numeric vector, \[0, 1\]: predicted probabilities
#'
#' @author EDG
#' @keywords internal
#' @noRd
brier_score <- function(true_int, predicted_prob) {
  true_int <- clean_int(true_int)
  check_float01inc(predicted_prob)
  mean((true_int - predicted_prob)^2)
} # /rtemis::brier_score

#' Convert labels to integers
#'
#' Convert factor labels to integers where the positive class is 1 and the negative class is 0.
#'
#' @param x Factor: True labels.
#' @param binclasspos Integer: Position of the factor level which is the positive class (binary classification only).
#'
#' @return Integer vector: 0, 1 where 1 is the positive class as defined by binclasspos.
#'
#' @author EDG
#' @keywords internal
#' @noRd
labels2int <- function(x, binclasspos = 2L) {
  stopifnot(is.factor(x))
  # Convert factor to 0, 1 where 1 is the positive class as defined by binclasspos
  if (binclasspos == 1L) {
    xi <- 2L - as.integer(x)
  } else {
    xi <- as.integer(x) - 1L
  }
  xi
} # /rtemis::labels2int

# classification_metrics() ----
#' Classification Metrics
#'
#' @details
#' Note that auc_method = "pROC" is the only one that will output an AUC even if
#' one or more predicted probabilities are NA.
#'
#' @param true_labels Factor: True labels.
#' @param predicted_labels Factor: predicted values.
#' @param predicted_prob Optional Numeric vector or matrix \[0, 1\]: Predicted
#'   probabilities, one column per class. A single column is the positive
#'   class's probability and enables the probability-based metrics (AUC, Brier
#'   score); a multiclass matrix has no single score per case, so those are not
#'   computed.
#' @param binclasspos Integer: Factor level position of the positive class in binary classification.
#' @param calc_auc Logical: If TRUE, calculate AUC. May be slow in very large datasets.
#' @param calc_brier Logical: If TRUE, calculate Brier score.
#' @param auc_method Character: "lightAUC", "pROC", "ROCR".
#' @param sample Optional Character \{"Training", "Validation", "Test",
#'   "Calibrated Training", "Calibrated Validation", "Calibrated Test"\}: Sample
#'   the metrics describe.
#' @param verbosity Integer: Verbosity level.
#'
#' @return `ClassificationMetrics` object.
#'
#' @author EDG
#' @export
#' @examples
#' # Assume positive class is "b"
#' true_labels <- factor(c("a", "a", "a", "b", "b", "b", "b", "b", "b", "b"))
#' predicted_labels <- factor(c("a", "b", "a", "b", "b", "a", "b", "b", "b", "a"))
#' predicted_prob <- c(0.3, 0.55, 0.45, 0.75, 0.57, 0.3, 0.8, 0.63, 0.62, 0.39)
#'
#' classification_metrics(true_labels, predicted_labels, predicted_prob)
#' classification_metrics(true_labels, predicted_labels, 1 - predicted_prob, binclasspos = 1L)
classification_metrics <- function(
  true_labels,
  predicted_labels,
  predicted_prob = NULL,
  binclasspos = 2L,
  calc_auc = TRUE,
  calc_brier = TRUE,
  auc_method = "lightAUC",
  sample = NULL,
  verbosity = 0L
) {
  # Checks ----
  # The probability-based metrics take one score per case, so a multiclass
  # matrix has nothing to give them and drops out here.
  predicted_prob <- positive_prob(predicted_prob)
  n_classes <- nlevels(true_labels)

  # Check same levels in
  if (!all(levels(true_labels) == levels(predicted_labels))) {
    rtemis.core::abort(
      "True and predicted labels must have the same levels, in the same order.",
      "\n     levels(true_labels): ",
      paste(levels(true_labels), collapse = ", "),
      "\nlevels(predicted_labels): ",
      paste(levels(predicted_labels), collapse = ", "),
      class = c("rtemis_level_mismatch", "rtemis_data_error")
    )
  }

  # Positive class ----
  # For confusion table, make positive class the first factor level
  if (n_classes == 2 && binclasspos == 2L) {
    true_labels <- factor(true_labels, levels = rev(levels(true_labels)))
    predicted_labels <- factor(
      predicted_labels,
      levels = rev(levels(predicted_labels))
    )
  }
  true_levels <- levels(true_labels)

  # Levels already set so that the first level is the positive class. NULL, not
  # NA, when there is no such thing: NULL is the only "unset" value, and a bare
  # NA would make the field logical for multiclass and character for binary.
  positive_class <- if (n_classes == 2L) true_levels[1L] else NULL
  if (verbosity > 0L) {
    if (n_classes == 2) {
      msg(
        "There are two outcome classes:",
        highlight(paste(rev(true_levels), collapse = ", "))
      )
      msg("        The positive class is:", highlight(positive_class))
    } else {
      msg(
        "There are",
        n_classes,
        "classes:",
        highlight(paste(rev(true_levels), collapse = ", "))
      )
    }
  }
  # Dimension names are not set here: `ClassificationMetrics` stores the long
  # form and derives the wide table, which labels both dimensions uniformly.
  tbl <- table(true_labels, predicted_labels)

  class <- list()
  overall <- list()
  class[["Totals"]] <- rowSums(tbl)
  class[["Predicted_totals"]] <- colSums(tbl)
  Total <- sum(tbl)
  class[["Hits"]] <- diag(tbl)
  # class[["Misses"]] <- class[["Totals"]] - class[["Hits"]]
  class[["sensitivity"]] <- class[["Hits"]] / class[["Totals"]]
  class[["Condition_negative"]] <- Total - class[["Totals"]]
  class[["True_negative"]] <- Total -
    class[["Predicted_totals"]] -
    (class[["Totals"]] - class[["Hits"]])
  class[["specificity"]] <- class[["True_negative"]] /
    class[["Condition_negative"]]
  class[["balanced_accuracy"]] <- .5 *
    (class[["sensitivity"]] + class[["specificity"]])
  # PPV = true positive / predicted condition positive
  class[["ppv"]] <- class[["Hits"]] / class[["Predicted_totals"]]
  # NPV  = true negative / predicted condition negative
  class[["npv"]] <- class[["True_negative"]] /
    (Total - class[["Predicted_totals"]])
  class[["f1"]] <- 2 *
    (class[["ppv"]] * class[["sensitivity"]]) /
    (class[["ppv"]] + class[["sensitivity"]])

  # Binary vs Multiclass ----
  if (n_classes == 2) {
    overall[["sensitivity"]] <- class[["sensitivity"]][1]
    overall[["specificity"]] <- class[["specificity"]][1]
    overall[["balanced_accuracy"]] <- class[["balanced_accuracy"]][1]
    overall[["ppv"]] <- class[["ppv"]][1]
    overall[["npv"]] <- class[["npv"]][1]
    overall[["f1"]] <- class[["f1"]][1]
  } else {
    overall[["balanced_accuracy"]] <- mean(class[["sensitivity"]])
    overall[["f1"]] <- mean(class[["f1"]])
  }
  overall[["accuracy"]] <- sum(class[["Hits"]]) / Total

  # Probability-based metrics ----
  if (!is.null(predicted_prob) && n_classes == 2L) {
    # Positive class has been set to first level
    true_int <- 2L - as.integer(true_labels)
    if (calc_auc) {
      overall[["auc"]] <- auc(
        true_int = true_int,
        predicted_prob = predicted_prob,
        method = auc_method
      )
    }
    if (calc_brier) {
      overall[["brier_score"]] <- brier_score(true_int, predicted_prob)
    }
    # overall[["Log loss"]] <- logloss(true_int, predicted_prob)
  }

  # Outro ----
  overall <- as.data.frame(do.call(cbind, overall))
  rownames(overall) <- NULL
  # `level` carries the outcome level each row describes. It is a column rather
  # than a row name because row names have no row-oriented JSON form: they would
  # be dropped on serialization, leaving unlabeled per-class metrics.
  class <- data.frame(
    level = rownames(tbl),
    sensitivity = class[["sensitivity"]],
    specificity = class[["specificity"]],
    balanced_accuracy = class[["balanced_accuracy"]],
    ppv = class[["ppv"]],
    npv = class[["npv"]],
    f1 = class[["f1"]],
    row.names = NULL
  )

  ClassificationMetrics(
    sample = sample,
    confusion_matrix = tbl,
    overall = overall,
    class = class,
    positive_class = positive_class
  )
} # /rtemis::classification_metrics


# regression_metrics() ----
#' Regression Metrics
#'
#' @param true Numeric vector: True values.
#' @param predicted Numeric vector: Predicted values.
#' @param na.rm Logical: If TRUE, remove NA values before computation.
#' @param sample Optional Character \{"Training", "Validation", "Test",
#'   "Calibrated Training", "Calibrated Validation", "Calibrated Test"\}: Sample
#'   the metrics describe.
#'
#' @return `RegressionMetrics` object.
#'
#' @author EDG
#' @export
#' @examples
#' true <- rnorm(100)
#' predicted <- true + rnorm(100, sd = 0.5)
#' regression_metrics(true, predicted)
regression_metrics <- function(
  true,
  predicted,
  na.rm = TRUE,
  sample = NULL
) {
  RegressionMetrics(
    mae = mae(true, predicted, na.rm = na.rm),
    mse = mse(true, predicted, na.rm = na.rm),
    rmse = rmse(true, predicted, na.rm = na.rm),
    rsq = rsq(true, predicted),
    sample = sample
  )
} # /rtemis::regression_metrics