rNIX/R/roc.R

#' Create a ROC curve
#'
#' @param labels vector of outcomes
#' @param scores vector of variable values
#' @param decreasing ordering of scores
#'
#' @return data frame with columns TPR, FPR, labels, scores, V 
#' V is combination of V0 (negative cases) and V1 (positive cases)
#' V0 is portion of positive cases with value larger than score
#' V1 is portion of negative cases with value below score
#'
#' @export

simple.roc <- function(labels, scores, decreasing = TRUE) {
  # Odstranimo NA vrednosti iz labels in scores
  valid_data = stats::complete.cases(labels, scores)
  labels = labels[valid_data]
  scores = scores[valid_data]

  labels <- labels[order(scores, decreasing=decreasing)]
  scores <- scores[order(scores, decreasing=decreasing)]

  labels = base::as.logical(labels)

  x = base::data.frame(TPR = base::cumsum(labels) / base::sum(labels),
                 FPR = base::cumsum(!labels) / base::sum(!labels),
                 labels, scores)

  s1 = scores[labels]
  s0 = scores[!labels]
  #V0 are 0 components of test, for each test-negative case find number of positive cases with value above it (ideally, count is equal to n1, number of positive cases)
  V0 = base::sapply(s0, psi01, s=s1)
  #V1 are 1 components of test; for each test-positive case find number of negative cases with score below its score (ideally, count is equal to n0, number of negative cases)
  V1 = base::sapply(s1, psi10, s=s0)
  n0=base::length(s0)
  n1=base::length(s1)
  #convert V0 to portions, to [0,1] range
  V0=V0/n1
  #convert V1 to portions, to [0,1] range
  V1=V1/n0

  V = labels
  V[labels] <- V1
  V[!labels] <- V0
  x$V = V
  x
}

psi01 <- function(y, s) {
  v1 = base::sum(s > y, na.rm = TRUE)
  v2 = base::sum(s == y, na.rm = TRUE)
  v1 + 0.5 * v2
}

psi10 <- function(x, s) {
  v1 = base::sum(s < x, na.rm = TRUE)
  v2 = base::sum(s == x, na.rm = TRUE)
  v1 + 0.5 * v2
}

#' Calculate AUC
#' @param roc as returned by simple.roc
#'
#' @return AUC (numeric)
#'
#' @export
simple.getAUC <- function(roc) {
  V0 = roc[!roc$labels, 'V']
  mean(V0)
}

#' Calculate metrics associated with ROC 
#'
#' @param ROC roc as returned by simple.roc
#'
#' @return a list with elemnts: : FPR, TPR, threshold, FPR CI95:min , FPR CI95: max, TPR CI95:min, TPR CI95:max
#' all evaluated at optimum point (Youden index)
#'
#' @export

simple.compute_roc_metrics <- function(ROC) {

#determine optimal point on ROC curve using Youden index
#ROC is the output of simple.roc
#output is a list : tpr, tpr CI95:min , tpr CI95: max, fpr, fpr CI95:min, fpr CI95:max, value of score,
#all evaluated at optimum point, 
   
   n=base::sum(ROC$labels)
   m=base::sum(!ROC$labels)
   nt=base::length(ROC$TPR)
   dist<-base::abs(ROC$TPR-ROC$FPR)
   imax<-base::which.max(dist)
   hTPR<-stats::prop.test(x=ROC$TPR[imax]*n,n=n,conf.level=0.95,correct=FALSE)
   hFPR<-stats::prop.test(x=ROC$FPR[imax]*m,n=m,conf.level=0.95,correct=FALSE)
   
  
   return(base::list(
    FPR = ROC$FPR[imax],
    TPR = ROC$TPR[imax],
    threshold = ROC$scores[imax],
    lFPR = hFPR$conf.int[1],
    hFPR = hFPR$conf.int[2],
    lTPR = hTPR$conf.int[1],
    hTPR = hTPR$conf.int[2]
  ))
}