R/cusTH_function.R

In snwikaij/GRASS: GRadient AnalySiS for ecologists

cusTH <- function(samp, tax, var, abu, conf.level = 0.95, nboot = 500){

  if(is.numeric(var) != T)
  {stop("var must be numeric.")}
  if(is.numeric(abu) != T)
  {stop("abu must be numeric.")}
  if(nboot < 100)
  {stop("nboot must be 100 or larger.")}
  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.")}

  df1        <- setNames(cbind.data.frame(samp, tax, var, abu),c("samp", "tax", "var", "abu"))
  df1        <- df1[df1$tax %in%  names(table(df1$tax))[table(df1$tax) >= 10],]

  custest    <- function(x, conf.level = conf.level, nboot = nboot, print = T){

    listb    <- list()
    for(b in 1:nboot){
      if(print == T){print(b)}
      df        <- x[sample(nrow(x), nrow(x), replace = T),]
      df        <- df[order(df$var),]

      tot       <- aggregate(data=df, tax~samp*var, length)
      tot$tax   <- 1
      tot$custax<- cumsum(tot$tax)/max(cumsum(tot$tax))

      totabu       <- aggregate(data=df, abu~samp*var, sum)
      totabu$cusabu<- cumsum(totabu$abu)/max(cumsum(totabu$abu))

      tot        <- cbind(tot, totabu)[-c(5,6)]
      tot$sumcus <- tot$custax+tot$cusabu
      tot$custot <- (tot$sumcus-min(tot$sumcus))/(max(tot$sumcus)-min(tot$sumcus))
      tot        <- tot[-c(7)]

      forloop   <- as.data.frame(matrix(ncol = 7, nrow = 0))

      for(i in unique(df$tax)){
        tax              <- df[df$tax == i,]
        tax$tax          <- 1

        taxsub           <- tot[tot$samp %in% tax$samp,]
        taxsub$taxsum    <- cumsum(taxsub$tax)/max(cumsum(taxsub$tax))
        taxsub$abusum    <- cumsum(taxsub$abu)/max(cumsum(taxsub$abu))
        taxsub$taxtot    <- taxsub$taxsum+taxsub$abusum
        taxsub$taxtot    <- (taxsub$taxtot-min(taxsub$taxtot))/(max(taxsub$taxtot)-min(taxsub$taxtot))
        rownames(taxsub) <- NULL

        loc              <- taxsub$var[which.max((taxsub$custot-taxsub$taxtot)^2)]
        resp             <- (taxsub$custot-taxsub$taxtot)[which.max((taxsub$custot-taxsub$taxtot)^2)]
        n                <- nrow(taxsub)
        auc              <- as.numeric(suppressWarnings(wilcox.test(taxsub$custot, taxsub$taxtot))[1])/n^2
        ming             <- as.numeric(quantile(taxsub$var, (1-conf.level)/2))
        maxg             <- as.numeric(quantile(taxsub$var, (1-conf.level)/2+conf.level))

        forloop[i,]      <- cbind(i, loc, ming, maxg, resp, auc, n)}
      listb[[b]]         <- forloop}

    dfb <- do.call(rbind.data.frame, listb)
    res <- setNames(cbind.data.frame(aggregate(data=dfb, V2~V1, function(x) mean(as.numeric(x), na.rm = T)),
                                     as.data.frame(aggregate(data=dfb, V2~V1, function(x) quantile(as.numeric(x), (1-conf.level)/2, na.rm = T)))[-1],
                                     as.data.frame(aggregate(data=dfb, V2~V1, function(x) quantile(as.numeric(x), (1-conf.level)/2+conf.level, na.rm = T)))[-1],
                                     aggregate(data=dfb, V3~V1, function(x) mean(as.numeric(x), na.rm = T))[-1],
                                     aggregate(data=dfb, V4~V1, function(x) mean(as.numeric(x), na.rm = T))[-1],
                                     aggregate(data=dfb, V5~V1, function(x) mean(as.numeric(x), na.rm = T))[-1],
                                     aggregate(data=dfb, V6~V1, function(x) mean(as.numeric(x), na.rm = T))[-1],
                                     aggregate(data=dfb, V6~V1, function(x) quantile(as.numeric(x), (1-conf.level)/2, na.rm = T))[-1],
                                     aggregate(data=dfb, V6~V1, function(x) quantile(as.numeric(x), (1-conf.level)/2+conf.level, na.rm = T))[-1],
                                     aggregate(data=dfb, V7~V1, function(x) mean(as.numeric(x), na.rm = T))[-1]),
                    c("Taxon", "Threshold", "LCI(t)", "HCI(t)", "min", "max", "Response",
                      "AUC", "LCI(a)", "HCI(a)", "n"))

    return(res)}

  results    <- custest(df1, conf.level = conf.level, nboot = nboot, print = T)}