R/AUCG.test_function.R

In snwikaij/GRASS: GRadient AnalySiS for ecologists

AUCG.test <- function(samp, tax, var, conf.level = 0.95, boot = FALSE, nboot = 500, pairwise = FALSE, mintax = 10){

  if(length(samp) != length(tax) ||
     length(tax)  != length(var) ||
     length(samp) != length(var))
  {stop("Vectors are not of the same length.")}
  if(is.numeric(var) != T)
  {stop("Vector var is not numeric.")}
  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.")}

  input             <- setNames(cbind.data.frame(samp, tax, var), c("samp", "tax", "var"))
  input             <- input[input$tax %in% names(table(input$tax))[table(input$tax) >= mintax],]
  input             <- na.omit(input)
  if(nrow(input) < 1)
  {stop("Not enough observations left after minsamp removal for analysis.")}

  quiet <- function(x){
    sink(tempfile())
    on.exit(sink())
    invisible(force(x))}

  dfsamp            <- input[!duplicated(input$samp),]
  outputlist        <- list()
  for(i in unique(input$tax)){
    dftax           <- input[input$tax == i,]
    dfnotax         <- dfsamp[!dfsamp$samp %in% c(as.character(dftax$samp)),]
    outputlist[[i]] <- quiet(AUC.test(dftax[,3], dfnotax[,3], conf.level = conf.level, boot = boot, nboot = nboot))}
  output            <- do.call(rbind.data.frame, outputlist)
  output            <- setNames(cbind(rownames(output), output), c("Taxon", colnames(output)))
  rownames(output)  <- NULL

  if(pairwise == TRUE){
    outputlist2         <- list()
    for(i in unique(input$tax)){
      dftax2            <- input[input$tax == i,]
      dfpair            <- as.data.frame(matrix(ncol=6, nrow=0))
      for(z in unique(na.omit(gsub(i,NA,input$tax)))){
        dftax3          <- input[input$tax == z,]
        dfpair[z,]      <- quiet(AUC.test(dftax2[,3], dftax3[,3], conf.level = conf.level))}
      outputlist2[[i]]  <- dfpair}

    output2             <- do.call(rbind.data.frame, outputlist2)
    outputnames         <- do.call(rbind, strsplit(as.character(rownames(output2)),".", fixed = T))
    output2             <- setNames(cbind(outputnames, output2),c("Taxon 1", "Taxon 2", colnames(output[-1])))
    rownames(output2)   <- NULL

    overlap.mat          <- reshape(output2[c("Taxon 1", "Taxon 2", "Overlap.estimate")], idvar="Taxon 1", timevar="Taxon 2", direction="wide")
    rownames(overlap.mat)<- overlap.mat$`Taxon 1`
    overlap.mat          <- overlap.mat[-1]
    overlap.mat          <- overlap.mat[c(ncol(my.mat),1:(ncol(overlap.mat)-1))]
    colnames(overlap.mat)<- rownames(overlap.mat)}

  Acom  <- sum(abs(output[2]-0.5)*2)/nrow(output)
  ntax  <- nrow(output)
  Acomr <- Acom/(sum(abs(cumsum(c(0, rep(1/(ntax-1),(ntax-1))))-0.5)*2)/ntax)

  comresp <- list(TAD    = Acom, TADr = Acomr, nr.tax = ntax,
                  `5%`   = as.numeric(quantile(var, 0.05, na.rm = T)),
                  `50%`  = as.numeric(quantile(var, 0.5, na.rm = T)),
                  `95%`  = as.numeric(quantile(var, 0.95, na.rm = T)),
                  `mean` = mean(var, na.rm = T))

  if(pairwise == TRUE){results <- list(comresp, output, output2, overlap.mat)
  names(results)<- c("Total assembly deviation", "Gradient deviation", "Pairwise deviation", "Overlap matrix")}
  else{results  <- list(comresp, output)
  names(results)<- c("Total assembly deviation", "Gradient deviation")}

  cat(paste0("TAD = ", round(results[[1]][[1]],2), "\n",
             "TADr = ", round(results[[1]][[2]],2), "\n",
             "n-taxa = ", ntax))
  invisible(results)}