R/performance_measures.R

In ipa-tys/ROCR: Visualizing the Performance of Scoring Classifiers

## ------------------------------------------------------------------------
## classical machine learning contingency table measures
## ------------------------------------------------------------------------

.performance.accuracy <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, (tn+tp) / length(predictions) )
  }

.performance.error.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, (fn+fp) / length(predictions) )
  }

.performance.false.positive.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, fp / n.neg )
  }

.performance.true.positive.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, tp / n.pos )
  }

.performance.false.negative.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, fn / n.pos )
  }

.performance.true.negative.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, tn / n.neg )
  }

.performance.positive.predictive.value <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    ppv <- tp / (fp + tp)
    list( cutoffs, ppv )
  }

.performance.negative.predictive.value <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    npv <- tn / (tn + fn)
    list( cutoffs, npv )
  }

.performance.prediction.conditioned.fallout <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    ppv <- .performance.positive.predictive.value(predictions, labels,
                                                  cutoffs, fp, tp, fn, tn,
                                                  n.pos, n.neg, n.pos.pred,
                                                  n.neg.pred)[[2]]
    list( cutoffs, 1 - ppv )
  }

.performance.prediction.conditioned.miss <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    npv <- .performance.negative.predictive.value(predictions, labels,
                                                  cutoffs, fp, tp, fn, tn,
                                                  n.pos, n.neg, n.pos.pred,
                                                  n.neg.pred)[[2]]
    list( cutoffs, 1 - npv )
  }

## ------------------------------------------------------------------------
## ...not actually performance measures, but very useful as a second axis
## against which to plot a "real" performance measure
## (popular example: lift charts)
## ------------------------------------------------------------------------

.performance.rate.of.positive.predictions <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, n.pos.pred / (n.pos + n.neg) )
  }

.performance.rate.of.negative.predictions <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list( cutoffs, n.neg.pred / (n.pos + n.neg) )
  }


## ------------------------------------------------------------------------
## Classical statistical contingency table measures
## ------------------------------------------------------------------------

.performance.phi <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    list(cutoffs,
         (tn*tp - fn*fp) /
           (sqrt(n.pos) * sqrt(n.neg) * sqrt(n.pos.pred) * sqrt(n.neg.pred))
        )
  }

.performance.mutual.information <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    n.samples <- n.pos + n.neg
    mi <- c()
    for (k in 1:length(cutoffs)) {
      kij <- rbind( c(tn[k],fn[k]), c(fp[k],tp[k]) )

      ki.j. <- rbind(c(n.neg * n.neg.pred[k], n.neg.pred[k] * n.pos),
                     c(n.neg * n.pos.pred[k], n.pos * n.pos.pred[k]))

      log.matrix <- log2( kij / ki.j.)
      log.matrix[kij/ki.j.==0] <- 0

      mi <- c(mi,  log2(n.samples) + sum( kij * log.matrix) / n.samples  )
    }

    list( cutoffs, mi )
  }

#' @importFrom stats chisq.test
.performance.chisq <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    chisq <- c()
    for (i in 1:length(cutoffs)) {
      A <- rbind( c( tn[i], fn[i]), c(fp[i], tp[i]) )
      chisq <- c(chisq, stats::chisq.test(A, correct=FALSE)$statistic )
    }
    list( cutoffs, chisq )
  }

.performance.odds.ratio <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {


    list( cutoffs, tp * tn / (fn * fp) )

  }

## ------------------------------------------------------------------------
## Other measures based on contingency tables
## ------------------------------------------------------------------------

.performance.lift <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    n.samples <- n.pos + n.neg
    list( cutoffs, (tp / n.pos) / (n.pos.pred / n.samples) )
  }

.performance.f <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, alpha) {

    prec <- .performance.positive.predictive.value(predictions, labels,
                                                   cutoffs, fp, tp, fn, tn,
                                                   n.pos, n.neg, n.pos.pred,
                                                   n.neg.pred)[[2]]
    list( cutoffs,  1/ ( alpha*(1/prec) + (1-alpha)*(1/(tp/n.pos))  ) )
  }

#' @importFrom grDevices chull
.performance.rocconvexhull <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    x <- fp / n.neg
    y <- tp / n.pos

    finite.bool <- is.finite(x) & is.finite(y)
    x <- x[ finite.bool ]
    y <- y[ finite.bool ]
    if (length(x) < 2) {
      stop("Not enough distinct predictions to compute ROC convex hull.")
    }

    ## keep only points on the convex hull
    ind <- grDevices::chull(x, y)
    x.ch <- x[ind]
    y.ch <- y[ind]

    ## keep only convex hull points above the diagonal, except (0,0)
    ## and (1,1)
    ind.upper.triangle <- x.ch < y.ch
    x.ch <- c(0, x.ch[ind.upper.triangle], 1)
    y.ch <- c(0, y.ch[ind.upper.triangle], 1)

    ## sort remaining points by ascending x value
    ind <- order(x.ch)
    x.ch <- x.ch[ind]
    y.ch <- y.ch[ind]

    list( x.ch, y.ch )
  }

## ----------------------------------------------------------------------------
## Cutoff-independent measures
## ----------------------------------------------------------------------------

#' @importFrom stats approxfun
.performance.auc <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, fpr.stop) {

    x <- fp / n.neg
    y <- tp / n.pos

    finite.bool <- is.finite(x) & is.finite(y)
    x <- x[ finite.bool ]
    y <- y[ finite.bool ]
    if (length(x) < 2) {
      stop(paste("Not enough distinct predictions to compute area",
                 "under the ROC curve."))
    }

    if (fpr.stop < 1) {
      ind <- max(which( x <= fpr.stop ))
      tpr.stop <- stats::approxfun( x[ind:(ind+1)], y[ind:(ind+1)] )(fpr.stop)
      x <- c(x[1:ind], fpr.stop)
      y <- c(y[1:ind], tpr.stop)
    }

    ans <- list()
    auc <- 0
    for (i in 2:length(x)) {
      auc <- auc + 0.5 * (x[i] - x[i-1]) * (y[i] + y[i-1])
    }
    ans <- list( c(), auc)
    names(ans) <- c("x.values","y.values")
    return(ans)
  }

# written by Thomas Unterthiner (unterthiner@bioinf.jku.at)
.performance.aucpr <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      
      
      tmp <- aggregate(list(fp=fp), by=list(tp=tp), min)
      tp <- tmp$tp
      fp <- tmp$fp
      prec <- tp / (fp + tp)
      rec <- tp / n.pos
      if (fp[1] == 0 & tp[1] == 0) {
          prec[1] = 1
      }

      finite.bool <- is.finite(prec) & is.finite(rec)
      prec <- prec[ finite.bool ]
      rec <- rec[ finite.bool ]	
      if (length(rec) < 2) {
          stop(paste("Not enough distinct predictions to compute area",
                     "under the Precision/Recall curve."))
      }
            
      # if two points are too distant from each other, we need to 
      # correctly interpolate between them. This is done according to
      # Davis & Goadrich, 
      #"The Relationship Between Precision-Recall and ROC Curves", ICML'06
      for (i in seq_along(rec[-length(rec)])) {
          if (tp[i+1] - tp[i] > 2) {
              skew = (fp[i+1]-fp[i]) / (tp[i+1]-tp[i])
              x = seq(1, tp[i+1]-tp[i], by=1)
              rec <- append(rec, (x+tp[i])/n.pos, after=i)
              prec <- append(prec, (x+tp[i])/(tp[i]+fp[i]+x+ skew*x), after=i)
          }
      }

      auc <- 0
      for (i in seq.int(from = 2, to = length(rec))) {
          auc <- auc + 0.5 * (rec[i] - rec[i-1]) * (prec[i] + prec[i-1])
      }
    
      ans <- list( c(), auc)
      names(ans) <- c("x.values","y.values")
      return(ans)
 }

#' @importFrom stats uniroot approxfun
.performance.precision.recall.break.even.point <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    pred <- prediction( predictions, labels)
    perf <- performance( pred, measure="prec", x.measure="rec")
    x <- rev(perf@x.values[[1]])
    y <- rev(perf@y.values[[1]])
    alpha <- rev(perf@alpha.values[[1]])

    finite.bool <- is.finite(alpha) & is.finite(x) & is.finite(y)
    x <- x[ finite.bool ]
    y <- y[ finite.bool ]
    alpha <- alpha[ finite.bool ]

    if (length(x) < 2) {
      stop(paste("Not enough distinct predictions to compute",
                 "precision/recall intersections."))
    }
    intersection.cutoff <- c()
    intersection.pr <- c()

    ## find all intersection points by looking at all intervals (i,i+1):
    ## if the difference function between x and y has different signs at the
    ## interval boundaries, then an intersection point is in the interval;
    ## compute as the root of the difference function
    if ( (x[1]-y[1]) == 0) {
      intersection.cutoff <- c( alpha[1] )
      intersection.pr <- c( x[1] )
    }

    for (i in (1:(length(alpha)-1))) {
      if ((x[i+1]-y[i+1]) == 0) {
        intersection.cutoff <- c( intersection.cutoff, alpha[i+1] )
        intersection.pr <- c( intersection.pr, x[i+1] )
      } else if ((x[i]-y[i])*(x[i+1]-y[i+1]) < 0 ) {
        ans <- stats::uniroot(stats::approxfun(c(alpha[i], alpha[i+1] ),
                                               c(x[i]-y[i], x[i+1]-y[i+1])),
                       c(alpha[i],alpha[i+1]))
        intersection.cutoff <- c(intersection.cutoff, ans$root)
        intersection.pr <- c(intersection.pr, ans$f.root)
      }
    }

    list( rev(intersection.cutoff), rev(intersection.pr) )
  }

#' @importFrom stats median
.performance.calibration.error <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, window.size) {

    if (window.size > length(predictions)) {
      stop("Window size exceeds number of predictions.")
    }
    if (min(predictions)<0 || max(predictions)>1) {
      stop("Calibration error needs predictions between 0 and 1")
    }

    pos.label <- levels(labels)[2]
    neg.label <- levels(labels)[1]

    ordering <- rev(order( predictions ))
    predictions <- predictions[ ordering ]
    labels <- labels[ ordering ]

    median.cutoffs <- c()
    calibration.errors <- c()

    for (left.index in 1 : (length(predictions) - window.size+1) ) {
      right.index <- left.index + window.size - 1
      pos.fraction <-
        sum(labels[left.index : right.index] == pos.label) / window.size
      mean.prediction <- mean( predictions[ left.index : right.index ] )

      calibration.errors <- c(calibration.errors,
                              abs(pos.fraction - mean.prediction))
      median.cutoffs <- c(median.cutoffs,
                          stats::median(predictions[left.index:right.index]))
    }
    list( median.cutoffs, calibration.errors )
  }


.performance.mean.cross.entropy <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    if (! all(levels(labels)==c(0,1)) ||
        any(predictions<0) || any(predictions>1) ) {
      stop(paste("Class labels need to be 0 and 1 and predictions between",
                 "0 and 1 for mean cross entropy."))
    }

    pos.label <- levels(labels)[2]
    neg.label <- levels(labels)[1]

    list( c(), - 1/length(predictions) *
            (sum( log( predictions[which(labels==pos.label)] ))  +
               sum( log( 1 - predictions[which(labels==neg.label)] ))) )
  }


.performance.root.mean.squared.error <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    ## convert labels from factor to numeric values
    labels <- as.numeric(levels(labels))[labels]
    if (any(is.na(labels))) {
      stop("For rmse predictions have to be numeric.")
    }
    list( c(),  sqrt( 1/length(predictions) *
                        sum( (predictions - labels)^2 ))  )
  }

## ----------------------------------------------------------------------------
## Derived measures:
## ----------------------------------------------------------------------------

.performance.sar <- function( predictions, labels, cutoffs, fp, tp, fn, tn,
                              n.pos, n.neg, n.pos.pred, n.neg.pred) {

  pred <- prediction( predictions, labels)
  perf.acc <- performance( pred, measure="acc")
  perf.rmse <- performance( pred, measure="rmse")
  perf.auc <- performance( pred, measure="auc")

  list(cutoffs,
       1/3 * (perf.acc@y.values[[1]] +
                (1 - perf.rmse@y.values[[1]]) +
                perf.auc@y.values[[1]]))
}

## ----------------------------------------------------------------------------
## Measures taking into account actual cost considerations
## ----------------------------------------------------------------------------

.performance.expected.cost <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    ## kick out suboptimal values (i.e. fpr/tpr pair for which another one
    ## with same fpr and higher tpr exists,
    ## or one for which one with same tpr but lower fpr exists

    if (n.neg==0 || n.pos==0) {
      stop(paste("At least one positive and one negative sample are",
                 "needed to compute a cost curve."))
    }
    fpr <- fp / n.neg
    tpr <- tp / n.pos

    ## sort by fpr (ascending), in case of ties by descending tpr
    ind <- order(fpr,-tpr)

    fpr <- fpr[ind]
    tpr <- tpr[ind]
    ## for tied fprs, only the one with the highest tpr is kept
    ind <- !duplicated(fpr)
    fpr <- fpr[ind]
    tpr <- tpr[ind]

    ## for tied tprs, only keep the one with the lowest fpr
    ind <- order(-tpr,fpr)
    fpr <- fpr[ind]
    tpr <- tpr[ind]
    ind <- !duplicated(tpr)
    fpr <- fpr[ind]
    tpr <- tpr[ind]

    if (!any(0==fpr & 0==tpr)) {
      fpr <- c(0,fpr)
      tpr <- c(0,tpr)
    }
    if (!any(1==fpr & 1==tpr)) {
      fpr <- c(fpr,1)
      tpr <- c(tpr,1)
    }

    ## compute all functions
    f <- list()
    for (i in 1:length(fpr)) {
      f <- c(f, .construct.linefunct( 0, fpr[i], 1, 1-tpr[i] ))
    }

    ## compute all intersection points
    x.values <- c()
    y.values <- c()
    for (i in 1:(length(fpr)-1)) {
      for (j in (i+1):length(fpr)) {
        ans <- .intersection.point( f[[i]], f[[j]] )
        if (all(is.finite(ans))) {
          y.values.at.current.x <- c()
          for (k in 1:length(f)) {
            y.values.at.current.x <- c(y.values.at.current.x,
                                       f[[k]](ans[1]))
          }
          if (abs(ans[2] - min(y.values.at.current.x )) <
              sqrt(.Machine$double.eps)) {

            x.values <- c(x.values, ans[1])
            y.values <- c(y.values, ans[2])
          }
        }
      }
    }

    if (!any(0==x.values & 0==y.values)) {
      x.values <- c(0,x.values)
      y.values <- c(0,y.values)
    }
    if (!any(1==x.values & 0==y.values)) {
      x.values <- c(x.values,1)
      y.values <- c(y.values,0)
    }

    ind <- order( x.values)
    list( x.values[ind], y.values[ind] )
  }


.performance.cost <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, cost.fp, cost.fn) {

    n.samples <- n.pos + n.neg
    cost <- ((n.pos / n.samples) * (fn / n.pos) * cost.fn +
               (n.neg / n.samples) * (fp / n.neg) * cost.fp)
    list( cutoffs, cost )
  }