R/test_fourth.R

In weimea: Weighted Mean Analysis

#' Testing the relationship species attributes and sample attributes by fourth corner method
#' 
#' Function calculating fourth corner analysis between sample attributes (environmental variables) and species attributes (traits) and testing it by row based, column based or max permutation test.
#' 
#' @name test_fourth
#' @param env Vector, matrix or data frame with sample attributes (environemntal variables).
#' @param com Matrix of data frame with species composition data.
#' @param traits Vector, matrix or data frame with species attributes (traits)
#' @param cwm An object of the class \code{cwm}.
#' @param perm Number of permutations.
#' @param test Vector of character values. Which test should be conducted? Partial match to \code{'rowbased'} row-based permutation test (called Model 2 in the Legendre et al. 1997 and Dray & Legendre 2008), \code{'colbased'} for column-based permutation test (Model 4), and \code{'max'} for max test (selecting the higher from \code{rowbased} and \code{colbased} result, Model 6). \code{'all'} includes all three tests. See \code{Details}.
#' @param adjustP Logical, default FALSE. Should be the P-values adjusted? If \code{adjustP = TRUE}, the last column in the results is adjusted by method selected in \code{p.adjust.method}.
#' @param p.adjust.method A string indicating the method of P-value adjustement, see \code{\link{p.adjust.methods}} for possible choices.
#' @param parallel NULL (default) or integer number. Number of cores for parallel calculation of modified permutation test. Maximum number of cores should correspond to number of available cores on the processor.
#' @param x,object object of the class \code{"cwm"}, generated by function \code{cwm}.
#' @param digits number of digits reported by \code{print} method on object of \code{"cwm"} class (default is 3).
#' @param missing.summary Logical; should be the summary of values missing in \code{env}, \code{cwm} and \code{traits} be printed along to the result of \code{test_cwm} analysis? Default is \code{TRUE}.
#' @param ... Other arguments for \code{print}, \code{summary} or \code{coef} functions (not implemented yet).

#' @export
#' @examples 
#' data (vltava)
#' vltava.env <- vltava$env[,c('pH', 'COVERE32')]
#' vltava.com <- vltava$herbs$spe
#' vltava.traits <- vltava$herbs$traits
#' fc <- test_fourth (env = vltava.env, com = vltava.com, traits = vltava.traits)
#' @details
#' Environmental variables, species composition and traits can be provided in two ways: either each object separately, using arguments \code{env}, \code{com} and \code{traits}, or using only \code{env} and \code{cwm}; in the later case, the function \code{extract} is used to extract the \code{com} and \code{traits} from the \code{cwm} object.
#' Specific issue related to weighted mean is the case of missing species attributes. In current implementation, species with missing species attributes are removed from sample x species matrix prior to permutation of species attributes among species. 
#' @return  Function \code{test_fourth} returns the list of the class \code{"testFOURTH"} (with \code{print} and \code{summary} methods), which contains the following items:
#' \itemize{
#'  \item \code{call} Function call.
#'  \item \code{out} Data frame with results of calculation; first two columns contain the original fourth corner metric (\code{r_fc}) and its Chessel transformation (\code{r_ch}); following columns contain some of the following results of permutation tests: \code{P.row}, \code{P.col} and \code{P.max}. If \code{adjustP = TRUE}, the last column contains P values adjusted by \code{p.adjust} function.
#'  \item \code{miss} List with missing values of \code{env}, \code{traits} or \code{cwm} (the last provided for compatibility with \code{test_cwm} function and refers to samples which, after removing missing \code{traits} values, do not contain any species.)
#'  \item \code{param} Setting of input parameters (\code{perm} and \code{test}).
#'  }
#' @seealso \code{\link{cwm}} \code{\link{test_cwm}}

#' @export
test_fourth <- function (env, com = NULL, traits = NULL, cwm = NULL, perm = 499, test = "max", parallel = NULL, p.adjust.method = p.adjust.methods, adjustP = FALSE)
{
  CALL <- match.call ()
  TEST <- c('none', 'rowbased', 'colbased', 'max', 'all')
  test <- match.arg (test, TEST, several.ok = T)
  env <- as.data.frame (as.matrix (env))
  if ( all (!is.null(cwm)&!is.null (com)&!is.null (traits)) | !is.null (cwm)&!is.null (com) | !is.null(cwm)&!is.null(traits)) stop ('Species composition and traits need to be provided using either com and traits argument or cwm argument; combining both is not possible!') 
  if (is.null (com) & is.null (traits)) {
    if (is.null (cwm)) stop ('You need to provide species composition data and traits, either using com and traits arguments, or cwm!') else
    {
      traits <- extract (cwm, 'traits')
      com <- extract (cwm, 'com')
    }
  }
  vars.names <- expand.grid (traits = colnames (traits), env = colnames (env), stringsAsFactors = FALSE)
  if (is.null (parallel)) {
      res <- apply (vars.names, 1, FUN = function (var12) test_fourthCpp (e = env[, var12['env'], drop = TRUE], L = as.matrix (com), t = traits[, var12['traits'], drop = TRUE], test = test, perm = perm))
    }
    else
    {
      cl <- parallel::makeCluster (parallel)
      parallel::clusterExport (cl, varlist = c('vars.names', 'com', 'traits', 'env', 'test', 'perm'), envir = environment ())
      parallel::clusterEvalQ (cl, eval (call ('library', 'weimea')))
      res <- parallel::parApply (cl, vars.names, 1, FUN = function (var12) test_fourthCpp (e = as.matrix (env[, var12['env'], drop = F]), L = as.matrix (com), t = as.matrix (traits[, var12['traits'], drop = F]), test = test, perm = perm))
    }
    names (res) <- apply (vars.names, 1, FUN = function (x) paste (x[1], x[2], sep = ' / '))

  res.coef <- lapply (res, FUN = function (x) x[1:(length (x)-3)])
  res.miss <- lapply (res, FUN = function (x) x[(length (x)-2):length (x)])
  result <- list (call = CALL, out = do.call (rbind.data.frame, res.coef), miss = res.miss, param = list (perm = perm, test = test))
  class (result) <- 'testFOURTH'
  return (result)
}

#' @rdname test_fourth
#' @export
print.testFOURTH <- function (x, digits = max(3, getOption("digits") - 3), missing.summary = FALSE, ...)
{
  cs.ind <- c(1, 2)
  tst.ind <- 0
  df.ind <- c(3, 4)
  cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), 
      "\n\n", sep = "")
  printCoefmat2 (x$out, digits = digits, cs.ind = cs.ind, tst.ind = tst.ind, df.ind = df.ind, ...)
  if (missing.summary) printMissum (x$miss)
  
}

printMissum <- function (x){
  miss <- lapply (x, FUN = function (a) lapply (a, length))
  miss <- as.matrix (do.call (rbind.data.frame, miss))
  colnames (miss) <- c('env', 'traits', 'samples')
  if (sum (miss) == 0) cat("\nThere are no missing env, traits or sample values.\n") else
  {
    cat("\nCount of missing values:\n")
    print.default (miss)
  }
  invisible (miss)
  }

#' @rdname test_fourth
#' @export
coef.testFOURTH <- function (object, ...)
  return (object$out)