R/AUC.test_function.R

In snwikaij/GRASS: GRadient AnalySiS for ecologists

AUC.test <- function(x, y, conf.level = 0.95, boot = FALSE, nboot = 500){

  if(is.numeric(x) != T)
  {stop("x must be numeric.")}
  if(is.numeric(y) != T)
  {stop("y must be numeric.")}
  if(length(x) == 0)
  {stop("Not enough samples for x, minimal 1.")}
  if(length(y) == 0)
  {stop("Not enough samples for y, minimal 1.")}
  if(conf.level > 1 || conf.level < 0 || is.numeric(conf.level) != T)
  {stop("Argument conf.level needs to be a numeric value inbetween 0 and 1.")}

  if(boot == F){
    W   <- suppressWarnings(stats::wilcox.test(x, y, conf.int = T, conf.level = conf.level))
    U   <- (length(x)*length(y)-W$statistic)
    A   <- (length(x)*length(y)-U)/(length(x)*length(y))
    O   <- 1-sqrt((as.numeric(A)-0.5)^2)*2

    Q1  <- A/(2-A)
    Q2  <- 2*A^2/(1+A)
    CI  <- sqrt((A*(1-A)+(length(x)-1)*(Q1-A^2)+(length(y)-1)*(Q2-A^2))/(length(x)*length(y)))
    LC  <- as.numeric(A-CI*stats::qnorm((1-conf.level)/2, lower.tail = F))
    if(LC < 0){LC <- 0}
    HC  <- as.numeric(A+CI*stats::qnorm((1-conf.level)/2, lower.tail = F))
    if(HC > 1){HC <- 1}
    res <- list(AUC.stat.estimate = as.numeric(A), Low.conf.AUC = LC, High.conf.AUC = HC,
                Overlap.estimate = as.numeric(O),
                x.obs = length(x), y.obs = length(y))
    rownames(res) <- NULL}else if(boot == T && nboot >= 100){
      strb <- as.numeric(NULL)
      for(i in 1:nboot){
        xb    <- sample(x, replace = T)
        yb    <- sample(y, replace = T)
        strb  <- c(strb, suppressWarnings(wilcox.test(xb,yb,conf.int = F)$stat[[1]]/(length(x)*length(y))))}

      A   <- mean(strb, na.rm = T)
      O   <- 1-sqrt((as.numeric(A)-0.5)^2)*2
      CI  <- quantile(strb, c((1-conf.level)/2, (1-conf.level)/2+conf.level), na.rm = T)
      LC  <- as.numeric(CI[1])
      if(LC < 0){LC <- 0}
      HC  <- as.numeric(CI[2])
      if(HC > 1){HC <- 1}
      res <- list(AUC.stat.estimate = A, Low.conf.AUC = LC, High.conf.AUC = HC,
                  Overlap.estimate = as.numeric(O),
                  x.obs = length(x), y.obs = length(y))}else{stop("Argument boot needs to be TRUE or FALSE and nboot >= 100.")}

  text <- paste(paste0("The area under the curve = ", round(A,2)), "\n",
                "At a confidence level of", paste0(conf.level*100,"%"), "\n",
                paste0("the lower confidence interval = ", round(LC,2)), "and", "\n",
                paste0("the higher confidence interval = ", round(HC,2)), "\n",
                paste0("The estimated sample overlap = ", round(O*100,2), "%" ))
  cat(text)
  invisible(res)}