R/boot.functions.R

In dkahle/TITAN2: Threshold Indicator Taxa Analysis

#' Calculate bootstrapped IndVal z scores
#'
#' This function implements resampling (with replacement) of the observed
#' environmental gradient and site-by-taxon matrix, and then calls the function
#' [getivz()] to obtain bootstrapped scores.
#'
#' Four pieces of information are obtained from every taxon during each
#' bootstrap replicate. If the argument 'imax' is TRUE, bootstrapped change
#' points are identified based on IndVal maxima, whereas if 'imax' is FALSE
#' z-score maxima are used instead.  In addition to the IndVal or z score
#' maxima, the value of the environmental gradient, the indicator direction, and
#' the p value are also retained for that point.
#'
#' In addition to the above metric matrix for each taxon (1), the z scores
#' across all candidate chnage points are retained from each replicate (2), as
#' well as the response direction (maxgrp) (3) and a sorted version of resampled
#' environmental values (4).  These four items are combined as a list object.
#'
#' @param bSeq An index used to determine the sequence number of the current
#'   bootstrap replicate.
#' @param env An environmental gradient.
#' @param taxa A site-by-taxon matrix of taxa counts at each sampling location.
#' @param ivTot A logical indicating whether IndVal scores should be calculated
#'   using total relative abundance or the mean relative abundace originally
#'   proposed by Dufrene and Legendre (1997). The default is to pass on the
#'   argument from the original TITAN funtion call.
#' @param minSplt The minimum bin size for partitioning along the environmental
#'   gradient.  The default is to pass on the argument from the original TITAN
#'   funtion call.
#' @param nPerm The number of replicates used by the permutation procedure (not
#'   to be confused with the number of bootstrap replicates). The default is to
#'   pass on the argument from the original TITAN funtion call.
#' @param memory A logical indicating whether scratch files should be used to
#'   store temporary data in order to preserve RAM during bootstrapping of large
#'   data sets.  The default is to pass on the argument from the original TITAN
#'   funtion call.
#' @param imax A logical indicating whether taxon-specific change points should
#'   be determined by IndVal maxima or z-score maxima (as in Baker and King
#'   2010). The default is to pass on the argument from the original TITAN
#'   funtion call.
#' @return A list of four elements: \describe{
#'
#'   \item{bt.metrics}{A matrix with nrow equal to number of taxa where the
#'   first column is the bootstrapped IndVal or z score maximum, the second is
#'   the environmental value, the third is the indicator direction, and the
#'   fourth is the p value at that point.}
#'
#'   \item{ivzs}{Z scores for all taxa across candidate change points in the
#'   replicate sample}
#'
#'   \item{bsrti}{A sorted version of the bootstrapped environmental gradient}
#'
#'   \item{rspdr}{Response direction (1 or 2) for all taxa across candidate
#'   change points in the replicate sample}
#'
#'   }
#' @references Baker, ME and RS King.  2010. A new method for detecting and
#'   interpreting biodiversity and ecological community thresholds. Methods in
#'   Ecology and Evolution 1(1): 25:37.
#' @seealso [getivz()], [ivzsums()]
#' @author M. Baker and R. King
#' @keywords TITAN bootstrap
tboot <- function(bSeq, env, taxa, ivTot = ivTot, minSplt = minSplt, nPerm = nPerm, memory = memory, imax = imax) {

  # below was replaced by the progress bar
  # cli_alert_info("{bSeq}.. ")

  # create empty arrays for storing sum(z-) and sum(z+) across envcls by replicate
  boot.metrics <- rep(NA, length(bSeq))
  numUnit <- length(env) # number along gradient
  rnum <- runif(1, 0, 100)
  max.btSumz <- rep(NA, 2)
  n_taxa <- ncol(taxa)

  for (i in 1:length(bSeq)) {

    # resample data with replacement, repeat part 1 preliminaries
    permuted_row_ndcs <- sample(1:numUnit, replace = TRUE) # generate permutation of row indices
    booted_taxa <- taxa[permuted_row_ndcs, ]
    permuted_env <- env[permuted_row_ndcs]
    bRankEnv <- rank(permuted_env, ties.method = "random")
    bSrti <- sort(permuted_env)

    bSrti2 <- sort(bRankEnv)
    bEnvCls <- bSrti[minSplt:(numUnit - minSplt)]
    nClass <- length(bEnvCls)
    boot.env <- matrix(NA, numUnit, nClass)
    for (c in 1:nClass) {
      boot.env[, c] <- (bRankEnv > bSrti2[minSplt + (c - 1)]) + 1
    }

    # create temporary matrix for storing z scores and group assignments
    ivzScores.bt <- getivz(boot.env, booted_taxa, ivTot, nPerm = nPerm, imax = imax, numClass = nClass)

    # replace NaNs with 0s
    for (k in (n_taxa + 1):(n_taxa * 3)) {
      for (j in 1:nClass) {
        if (is.nan(ivzScores.bt[k, j]))  ivzScores.bt[k, j] <- 0
      }
    }


    bt.metrics <- matrix(NA, n_taxa, 4)  #output matrix of metrics

    # obtain env level for IndVal maxima and group assignment for
    # each replicate bt.metric 1:4 include maxgrp, env.cp (either
    # imax or zmax), z.score, iv.p
    for (j in 1:n_taxa) {
      if (imax == FALSE) {
        maxSplt <- which.max(abs(ivzScores.bt[j + (n_taxa),]))
      } else {
        maxSplt <- which.max(abs(ivzScores.bt[j + (n_taxa * 2), ]))
      }
      bt.metrics[j, 1] <- ivzScores.bt[j, maxSplt]
      bt.metrics[j, 2] <- (bSrti[minSplt + maxSplt - 1] + bSrti[minSplt + maxSplt])/2
      bt.metrics[j, 3] <- ivzScores.bt[j + n_taxa, maxSplt]
      bt.metrics[j, 4] <- ivzScores.bt[(j + (n_taxa * 3)), maxSplt]
    }

    # obtain z scores and response direction
    ivzs <- ivzScores.bt[(n_taxa+1):(2*n_taxa),]
    rspdr <- ivzScores.bt[1:n_taxa,]
    if (memory) {
      if (!file.exists("temp.dir")) dir.create("temp.dir")
      write.table(ivzs, glue("temp.dir/boot.z.{rnum}.{i}.txt"), append = F, row.names = F, col.names = F)
      write.table(rspdr, glue("temp.dir/boot.rd.{rnum}.{i}.txt"), append = F, row.names = F, col.names = F)
      write.table(bSrti, glue("temp.dir/bSrti.{rnum}.{i}.txt"), append = F, row.names = F, col.names = F)
      bt.list <- list(bt.metrics, rnum)
    } else {
      bt.list <- list(bt.metrics, ivzs, bSrti, rspdr)
    }
    boot.metrics[i] <- list(bt.list)
  }

  boot.metrics
}






























#' Controls the allocation of bootstrap replicates
#'
#' A wrapper function for controlling the implementation of bootstrap replicates
#' using the function 'tboot' by sequential or parallel processing.
#'
#' If 'ncpus'>1 evaluates to TRUE, the function employs the package 'snow' to
#' implement parallel processing on multicore processors common in modern
#' desktop computers.  With some minor modification it is possible to configure
#' this code to allocate processes to cores on a high-performance computing
#' cluster (i.e., a supercomputer).  If 'ncpus'>1 evaluates to FALSE, the
#' function uses 'lapply' to run 'tboot' in sequence 1:nBoot times.
#'
#' @param env A vector of values for each sampling location along the
#'   environmental gradient.
#' @param taxa A site-by-taxon matrix of taxa counts at each sampling location.
#' @param ivTot A logical indicating whether IndVal scores should be calculated
#'   using total relative abundance or the mean relative abundace originally
#'   proposed by Dufrene and Legendre (1997). The default is to pass on the
#'   argument from the original TITAN funtion call.
#' @param boot A logical specifying whether or not to implement TITAN's'
#'   boostrap procedure. The default is to use the value specified in the
#'   original TITAN function call.
#' @param ncpus An argument specifying the number of processing cores used by
#'   the TITAN function call.  If ncpus>1 then parallel processing is
#'   implemented.  The default is to use the value specified in the original
#'   TITAN function call.
#' @param nBoot An argument specifying the number of bootstrap replicates.  The
#'   default is to use the value specified in the original TITAN function call.
#' @param minSplt An argument specifying minimum split size for partitioning
#'   along the environmental gradient.  The default is to use the value
#'   specified in the original TITAN function call.
#' @param nPerm The number of replicates used by the permutation procedure (not
#'   to be confused with the number of bootstrap replicates).
#' @param memory A logical indicating whether scratch files should be used to
#'   store temporary data in order to conserve active memory during
#'   bootstrapping of large data sets.  The default is to pass on the argument
#'   from the original TITAN funtion call.
#' @param imax A logical indicating whether taxon-specific change points should
#'   be determined by IndVal maxima or z-score maxima (as in TITAN v1.0). The
#'   default is to pass on the argument from the original TITAN funtion call.
#' @param numUnit An argument specifying the number of values along the
#'   environmental gradient.
#' @return A list of two items: \describe{
#'
#'   \item{bSeq}{An index of the sequence of bootstrap replicates. The structure
#'   of bSeq will differ for sequential or parallel processing.}
#'
#'   \item{ivz.bt.list}{Itself a list of four items comprising output passed on
#'   from function [tboot()]}
#'
#'   }
#' @references Baker, ME and RS King.  2010. A new method for detecting and
#'   interpreting biodiversity and ecological community thresholds. Methods in
#'   Ecology and Evolution 1(1): 25:37.
#' @seealso [tboot()], [small.boot()], [big.boot()], [titan()]
#' @author M. Baker and R. King
#' @keywords TITAN bootstrap
boot.titan <- function(env, taxa, ivTot = ivTot, boot = boot, ncpus = ncpus,
  nBoot = nBoot, minSplt = minSplt, nPerm = 250, memory = memory,
  imax = imax, numUnit = numUnit) {

  # if multiple cores are available, take advantage of parallel processing
  if (ncpus > 1) {
    cli_alert_info("Bootstrap resampling in parallel using {ncpus} CPUs... no index will be printed to screen.")
    # if (ncpus > nBoot) {
    #   cli_alert_info("Decreasing number of CPUs to number of bootstrap replicates ({nBoot}).")
    #   ncpus <- nBoot
    # }
    # cl <- parallel::makeCluster(rep("localhost", ncpus), type = "SOCK")
    cl <- parallel::makeCluster(rep("localhost", ncpus))
    bSeq <- parallel::clusterSplit(cl, 1:nBoot)
    ivz.bt.list <- parallel::clusterApply(cl, bSeq, tboot, env = env, taxa = taxa, ivTot = ivTot, minSplt = minSplt, nPerm = nPerm, memory = memory, imax = imax)
    parallel::stopCluster(cl)
  } else {
    # otherwise run bootstrap in sequence
    cli_alert_info("Bootstrap resampling in sequence...")
    bSeq <- 0
    # ivz.bt.list = lapply(1:nBoot, tboot, env = env, taxa = taxa, ivTot = ivTot, minSplt = minSplt, nPerm = nPerm, memory = memory, imax = imax)
    cli_progress_bar("Resampling", total = nBoot)
    ivz.bt.list <- vector(mode = "list", length = nBoot)
    for (i in 1:nBoot) {
      ivz.bt.list[[i]] <- tboot(i, env = env, taxa = taxa, ivTot = ivTot, minSplt = minSplt, nPerm = nPerm, memory = memory, imax = imax)
      cli_progress_update()
    }
    cli_progress_done()
  }

  list(bSeq, ivz.bt.list)
}























#' Summarizes raw output from TITAN's bootstrap procedure
#'
#' A function to take output from TITAN's bootstrap procedure and process it for
#' summary output.  The default is to perform this processing entirely within
#' active memory, but in the event of overflowing system capacity, an optional
#' program writes temporary files to a scratch directory to circumvent memory
#' limits.
#'
#' Use of 'small.boot' versus 'big.boot' is controlled by the argument 'memory'
#' in the original TITAN function call and passed to the wrapper function
#' 'titan'.  The two progams have identical functionality, but they accomplish
#' those functions differently to deal with memory limitations.
#'
#' For sequential processing of the bootsrtap, the index 'bSeq' is simply a
#' sequence from 1:nBoot that is printed to the screen. For parallel processing,
#' 'bSeq' is a list of length equal to 'ncpus', where each item is a segment of
#' the sequence allocated to each processing core.  Thus, depending on whether
#' 'ncpus'>1, the value of 'bSeq' is used differently to extract values from the
#' bootstrap output list.
#'
#' The first part of each function consists of defining output matrices, the
#' second involves extraction of output from the bootstrap list, the third part
#' involves calculating purity, reliability, the median z score, and quantiles
#' of the bootstrapped change points for each taxon.  These values are used to
#' complete the 'sppmax' output table and to identify the taxa that meet purity
#' and reliability criteria.  The final portion of each function finds the
#' maximum sum(z-), sum(z+), f.sum(z-), and f.sum(z+) for each bootstrap
#' replicate for later estimation of confidence intervals.  The final portion of
#' the summary involves calculating the filtered and unfiltered sum(z) scores
#' for each bootstrap replicate from the matrix of z scores and response
#' directions passed from the function boot.titan() within ivz.bt.list
#'
#' @param ivz.bt.list A list of output from each bootstrap replicate passed from
#'   [boot.titan()].
#' @param bSeq An index of the sequence of bootstrap replicates.
#' @param sppmax A taxon-specific summary output table for TITAN.
#' @param obs1 A binary vector indicating membership in the decreasing group of
#'   taxa.
#' @param obs2 A binary vector indicating membership in the increasing group of
#'   taxa.
#' @param nBoot An argument specifying the number of bootstrap replicates.  The
#'   default is to use the value specified in the original TITAN function call.
#' @param numClass An argument specifying the number of candidate partitions
#'   along the environmental gradient.
#' @param numUnit An argument specifying the number of values along the
#'   environmental gradient.
#' @param ncpus An argument specifying the number of processing cores used by
#'   the TITAN function call.  If ncpus>1 then parallel processing is indicated.
#'   The default is to use the value specified in the original TITAN function
#'   call.
#' @param pur.cut An argument specifying the cutoff value for determining
#'   purity.  The default is to use the value specified in the original TITAN
#'   function call.
#' @param rel.cut An argument specifying the cutoff value for determining
#'   reliability.  The default is to use the value specified in the original
#'   TITAN function call.
#' @param minSplt An argument specifying minimum split size of partitioning
#'   along the environmental gradient.  The default is to use the value
#'   specified in the original TITAN function call.
#' @return A list of six items: \describe{
#'
#'   \item{sppSub1}{A vector of taxon index numbers for pure and reliable
#'   decreasers}
#'
#'   \item{sppSub2}{A vector of taxon index numbers for pure and reliable
#'   increasers}
#'
#'   \item{sppmax}{The completed taxon-specific summary output table for TITAN}
#'
#'   \item{maxSumz}{A 2-column matrix of environmental values at sum(z-) and
#'   sum(z+) maxima across all bootstrap replicates}
#'
#'   \item{maxFsumz}{A 2-column matrix of environmental values at filtered
#'   sum(z-) and sum(z+) maxima across all bootstrap replicates}
#'
#'   \item{metricArray}{An array of group membership, env change points, z
#'   scores, and p values for passing to 'plot.IVecdf'}
#'
#'   }
#' @references Baker, ME and RS King.  2010. A new method for detecting and
#'   interpreting biodiversity and ecological community thresholds. Methods in
#'   Ecology and Evolution 1(1): 25:37.
#' @references Baker ME and RS King. 2013. Of TITAN and straw men: an appeal for
#'   greater understanding of community data. Freshwater Science 32(2):489-506.
#' @seealso [boot.titan()], [tboot()], [titan()]
#' @author M. Baker and R. King
#' @keywords TITAN purity reliability sum(z)
#' @name smallBigBoot






#' @rdname smallBigBoot
small.boot <- function(ivz.bt.list, bSeq, sppmax, obs1, obs2, nBoot,
  numClass, numUnit, ncpus, pur.cut, rel.cut, minSplt) {

  nTaxa <- nrow(sppmax)

  # create empty matrix for storing IndVal maxima and group assignments
  metricArray <- array(NA, c(nTaxa, 4, nBoot))

  # create empty matrices for storing nuids and z values at IndVal maxima for bootreps
  zArray     <- array(NA, c(nTaxa, numClass, nBoot))
  bEnvMatrix <- array(NA, c(nBoot, numUnit))
  rspdir     <- array(NA, c(nTaxa, numClass, nBoot))

  # create empty matrices to store env values of filtered and unfiltered sum(z) maxima by bootstrap replicate
  sumzBoot <- sumzBoot.f <- maxSumz <- maxFsumz <- matrix(NA, nBoot, 2)

  # create an array across all bootreps from parallel or sequenced result list
  # aseq <- 0
  if (ncpus > 1) {

    # combine lists of output into one list (from ncpus list of lists)
    ivz.bt.list <- do.call(c, ivz.bt.list)[unlist(bSeq)]
    for (k in seq_along(ivz.bt.list)) {
      metricArray[, , k] <- unlist(ivz.bt.list[[k]][[1]])
           zArray[, , k] <- unlist(ivz.bt.list[[k]][[2]])
         bEnvMatrix[k, ] <- unlist(ivz.bt.list[[k]][[3]])
           rspdir[, , k] <- unlist(ivz.bt.list[[k]][[4]])
    }

    # for (s in 1:ncpus) {
    #   # for (l in 1:length(bSeq[[s]])) {
    #   for (l in seq_along(bSeq[[s]])) {
    #     aseq <- aseq + 1
    #     metricArray[, , aseq] <- unlist(ivz.bt.list[[s]][[l]][[1]])
    #          zArray[, , aseq] <- unlist(ivz.bt.list[[s]][[l]][[2]])
    #      bEnvMatrix[aseq, ]   <- unlist(ivz.bt.list[[s]][[l]][[3]])
    #          rspdir[, , aseq] <- unlist(ivz.bt.list[[s]][[l]][[4]])
    #   }
    # }

  } else {

    for (i in 1:nBoot) {
      metricArray[, , i] <- unlist(ivz.bt.list[[i]][[1]][[1]])
           zArray[, , i] <- unlist(ivz.bt.list[[i]][[1]][[2]])
         bEnvMatrix[i, ] <- unlist(ivz.bt.list[[i]][[1]][[3]])
           rspdir[, , i] <- unlist(ivz.bt.list[[i]][[1]][[4]])
    }

  }

  purity1  <- rowMeans(metricArray[, 1, , drop=FALSE] == 1, na.rm = T)
  purity2  <- rowMeans(metricArray[, 1, , drop=FALSE] == 2, na.rm = T)
  reliab   <- rowMeans(metricArray[, 4, , drop=FALSE] < 0.05, na.rm = T)
  z.median <-    apply(metricArray[, 3, , drop=FALSE], 1, median, na.rm = T)
  quantiles <- t(apply(metricArray[, 2, , drop=FALSE], 1, quantile, probs = c(0.05, 0.1, 0.5, 0.9, 0.95), na.rm = T))
  sppmax[, 14] <- reliab
  sppmax[, 15] <- z.median
  sppmax[which(purity1 >= pur.cut & reliab >= rel.cut & obs1), 16] <- 1
  sppmax[which(purity2 >= pur.cut & reliab >= rel.cut & obs2), 16] <- sppmax[which(purity2 >= pur.cut & reliab >= rel.cut & obs2), 16] + 2
  sppmax[which(sppmax[, 4] == 1), 13] <- purity1[which(sppmax[, 4] == 1)]
  sppmax[which(sppmax[, 4] == 2), 13] <- purity2[which(sppmax[, 4] == 2)]
  sppmax[, 8:12] <- quantiles
  sppSub1 <- which(sppmax[, 16] == 1)
  sppSub2 <- which(sppmax[, 16] == 2)

  for (i in 1:nBoot) {
    sumzBoot[i, 1] <- which.max(colSums(zArray[, , i] * (rspdir[, , i] == 1), na.rm = T))
    if (length(sppSub1) > 1) {
      sumzBoot.f[i, 1] <- which.max(colSums(zArray[sppSub1, , i] * (rspdir[sppSub1, , i] == 1), na.rm = T))
    } else {
      if (length(sppSub1) > 0) {
        sumzBoot.f[i, 1] <- which.max(zArray[sppSub1, , i] * (rspdir[sppSub1, , i] == 1))
      }
    }
    sumzBoot[i, 2] <- which.max(colSums(zArray[, , i] * (rspdir[, , i] == 2), na.rm = T))
    if (length(sppSub2) > 1) {
      sumzBoot.f[i, 2] <- which.max(colSums(zArray[sppSub2, , i] * (rspdir[sppSub2, , i] == 2), na.rm = T))
    } else {
      if (length(sppSub2) > 0) {
        sumzBoot.f[i, 2] <- which.max(zArray[sppSub2, , i] * (rspdir[sppSub2, , i] == 2))
      }
    }
     maxSumz[i, 1] <- (bEnvMatrix[i,   (minSplt + sumzBoot[i, 1]) - 1] + bEnvMatrix[i,   (minSplt + sumzBoot[i, 1])])/2
     maxSumz[i, 2] <- (bEnvMatrix[i,   (minSplt + sumzBoot[i, 2]) - 1] + bEnvMatrix[i,   (minSplt + sumzBoot[i, 2])])/2
    maxFsumz[i, 1] <- (bEnvMatrix[i, (minSplt + sumzBoot.f[i, 1]) - 1] + bEnvMatrix[i, (minSplt + sumzBoot.f[i, 1])])/2
    maxFsumz[i, 2] <- (bEnvMatrix[i, (minSplt + sumzBoot.f[i, 2]) - 1] + bEnvMatrix[i, (minSplt + sumzBoot.f[i, 2])])/2
  }
  list(sppSub1, sppSub2, sppmax, maxSumz, maxFsumz, metricArray)
}















#' @rdname smallBigBoot
big.boot <- function(ivz.bt.list, bSeq, sppmax, obs1, obs2, nBoot, numClass, numUnit, ncpus, pur.cut, rel.cut, minSplt) {

  numTxa <- nrow(sppmax)
  # create empty matrix for storing IndVal maxima and group assignments
  metricArray <- array(NA, c(numTxa, 4, nBoot))

  # store env values of filtered and unfiltered sum(z) maxima by bootstrap replicate
  sumzBoot <- matrix(NA, nBoot, 2)
  sumzBoot.f <- matrix(NA, nBoot, 2)
  maxSumz <- matrix(NA, nBoot, 2)
  maxFsumz <- matrix(NA, nBoot, 2)
  track.seq <- rep(NA, nBoot)
  aseq = 0

  # create an array across all bootreps from parallel or sequenced result list
  if (ncpus > 1) {
    for (s in 1:ncpus) {
      for (l in 1:length(bSeq[[s]])) {
        aseq = aseq + 1
        metricArray[, , aseq] <- unlist(ivz.bt.list[[s]][[l]][[1]])
        track.seq[aseq] <- ivz.bt.list[[s]][[l]][[2]]
      }
    }
  } else {
    for (i in 1:nBoot) {
      aseq = aseq + 1
      metricArray[, , i] <- unlist(ivz.bt.list[[i]][[1]][[1]])
      track.seq[aseq] <- ivz.bt.list[[i]][[1]][[2]]
    }
  }

  purity1 <- rowMeans(metricArray[,1,] == 1, na.rm = T)
  purity2 <- rowMeans(metricArray[,1,] == 2, na.rm = T)
  reliab  <- rowMeans(metricArray[,4,] < 0.05, na.rm = T)
  z.median <- apply(metricArray[, 3, ], 1, median, na.rm = T)
  quantiles <- t(apply(metricArray[, 2, ], 1, quantile, probs = c(0.05, 0.1, 0.5, 0.9, 0.95), na.rm = T))

  sppmax[, 14] <- reliab
  sppmax[, 15] <- z.median
  sppmax[which(purity1 >= pur.cut & reliab >= rel.cut & obs1), 16] <- 1
  sppmax[which(purity2 >= pur.cut & reliab >= rel.cut & obs2), 16] <- sppmax[which(purity2 >= pur.cut & reliab >= rel.cut & obs2), 16] + 2
  sppmax[which(sppmax[, 4] == 1), 13] <- purity1[which(sppmax[,4] == 1)]
  sppmax[which(sppmax[, 4] == 2), 13] <- purity2[which(sppmax[,4] == 2)]
  sppmax[, 8:12] <- quantiles
  sppSub1 <- which(sppmax[, 16] == 1)
  sppSub2 <- which(sppmax[, 16] == 2)

  # read z score matrices in one at a time and obtain filtered sumz maxima
  if (ncpus > 1) {
    aseq = 0
    for (s in 1:ncpus) {
      for (l in 1:length(bSeq[[s]])) {
        aseq = aseq + 1
        ztab    = read.table(paste("temp.dir/boot.z.", track.seq[aseq], ".", l, ".txt", sep = ""))
        rspdir  = read.table(paste("temp.dir/boot.rd.", track.seq[aseq], ".", l, ".txt", sep = ""))
        bSrtEnv = read.table(paste("temp.dir/bSrti.", track.seq[aseq], ".", l, ".txt", sep = ""))
        sumzBoot[aseq, 1] <- which.max(colSums(ztab * (rspdir == 1), na.rm = T))

        if (length(sppSub1) > 1) {
          sumzBoot.f[aseq, 1] <- which.max(colSums(ztab[sppSub1,] * (rspdir[sppSub1, ] == 1), na.rm = T))
        } else {
          sumzBoot.f[aseq, 1] <- which.max(ztab[sppSub1,] * (rspdir[sppSub1, ] == 1))
        }
        sumzBoot[aseq, 2] <- which.max(colSums(ztab * (rspdir == 2), na.rm = T))
        if (length(sppSub2) > 1) {
          sumzBoot.f[aseq, 2] <- which.max(colSums(ztab[sppSub2,] * (rspdir[sppSub2, ] == 2), na.rm = T))
        } else {
          sumzBoot.f[aseq, 2] <- which.max(ztab[sppSub2,] * (rspdir[sppSub2, ] == 2))
        }

        maxSumz[aseq, 1]  <- (bSrtEnv[(minSplt +   sumzBoot[aseq,1]) - 1, ] + bSrtEnv[(minSplt +   sumzBoot[aseq,1]), ])/2
        maxSumz[aseq, 2]  <- (bSrtEnv[(minSplt +   sumzBoot[aseq,2]) - 1, ] + bSrtEnv[(minSplt +   sumzBoot[aseq,2]), ])/2
        maxFsumz[aseq, 1] <- (bSrtEnv[(minSplt + sumzBoot.f[aseq,1]) - 1, ] + bSrtEnv[(minSplt + sumzBoot.f[aseq,1]), ])/2
        maxFsumz[aseq, 2] <- (bSrtEnv[(minSplt + sumzBoot.f[aseq,2]) - 1, ] + bSrtEnv[(minSplt + sumzBoot.f[aseq,2]), ])/2
      }
    }
  } else {
    for (i in 1:nBoot) {
      ztab = read.table(paste("temp.dir/boot.z.", track.seq[i], ".", 1, ".txt", sep = ""))
      rspdir = read.table(paste("temp.dir/boot.rd.", track.seq[i], ".", 1, ".txt", sep = ""))
      bSrtEnv = read.table(paste("temp.dir/bSrti.", track.seq[i], ".", 1, ".txt", sep = ""))
      sumzBoot[i, 1] <- which.max(colSums(ztab * (rspdir == 1), na.rm = T))
      sumzBoot[i, 2] <- which.max(colSums(ztab * (rspdir == 2), na.rm = T))
      sumzBoot.f[i, 1] <- which.max(colSums(ztab[sppSub1,,i] * (rspdir[sppSub1,,i] == 1), na.rm = T))
      sumzBoot.f[i, 2] <- which.max(colSums(ztab[sppSub2,,i] * (rspdir[sppSub2,,i] == 2), na.rm = T))
       maxSumz[i, 1] <- (bSrtEnv[(minSplt +   sumzBoot[i,1]) - 1,] + bSrtEnv[(minSplt +   sumzBoot[i,1]),])/2
       maxSumz[i, 2] <- (bSrtEnv[(minSplt +   sumzBoot[i,2]) - 1,] + bSrtEnv[(minSplt +   sumzBoot[i,2]),])/2
      maxFsumz[i, 1] <- (bSrtEnv[(minSplt + sumzBoot.f[i,1]) - 1,] + bSrtEnv[(minSplt + sumzBoot.f[i,1]),])/2
      maxFsumz[i, 2] <- (bSrtEnv[(minSplt + sumzBoot.f[i,2]) - 1,] + bSrtEnv[(minSplt + sumzBoot.f[i,2]),])/2
    }
  }
  list(sppSub1, sppSub2, sppmax, maxSumz, maxFsumz, metricArray)
}