Indicator species (indicator value) functions

indvals.R

# Indicator value functions -----------------------------------------------

#' Returns the indicator value (Dufrene, 1997) for a given row of species counts
#' along with a vector of class assignments.
#' @param row:     vector of counts (usually a row in a counts matrix)
#' @param class:   which level in the grouping variable to test
#' @param classes: factor describing the grouping of the counts vector
.indval <- function(row, class, classes) {
  idxs <- classes == class
  A.ij <- sum(row[idxs]) / sum(row)
  B.ij <- sum(row[idxs] > 0) / sum(idxs)
  return(A.ij*B.ij)
}

#' Returns the indicator value and the p-value (as determined by permutation) 
#' for each feature (row) in a species matrix.
#' @param mat:   a species matrix (features=rows, samples=columns)
#' @param class: which level in the grouping vector to test
#' @param classes: a factor describing the grouping of the samples in the matrix
#' @param n.permutations: number of permutations to do to calculate the p value
#' @param p.adjust.method: which method to use to account for multiple testing
#'   (default: 'fdr')
get_indvals <- function(cts, group.factor, test.group, n.permutations=10000, p.adjust.method='fdr') {
  message("Performing indicator species permutation test")
  stopifnot(test.group %in% levels(group.factor))
  stopifnot(ncol(cts) == length(group.factor))
  cl <- makePSOCKcluster(detectCores())
  indvals <- adply(cts, 1, function(x) {
    indval <- .indval(x, test.group, group.factor)
    rand.indvals <- parSapply(cl=cl, c(1:n.permutations), function(y, test.group, group.factor, .indval){
      .indval(x, test.group, sample(group.factor))
    }, test.group, group.factor, .indval)
    p.value <- sum(rand.indvals >= indval) / length(rand.indvals)
    return(data.frame(indval=indval, p.value=p.value))
  })
  stopCluster(cl)
  indvals$adj.p.value <- p.adjust(indvals$p.value, method=p.adjust.method)
  indvals$otu <- indvals$X1
  indvals$X1 <- NULL
  return(indvals)
}