R/anova_sites.R

In douconca: Double Constrained Correspondence Analysis for Trait-Environment Analysis in Ecology

#' @title Utility function: community-level permutation test in Double 
#' Constrained Correspondence Analysis (dc-CA)
#'
#' @description
#' \code{anova_sites} performs the community-level permutation test of dc-CA 
#' when site weights vary.
#' The test uses residual predictor permutation (ter Braak 2022), which is 
#' robust against differences in site total abundance in the \code{response} 
#' in \code{\link{dc_CA}} (ter Braak & te Beest, 2022).
#' The arguments of the function are similar to those of 
#' \code{\link[vegan]{anova.cca}}, but more restricted. With equal site-totals 
#' as in \code{\link{dc_CA}}, \code{anova(object$RDAonEnv)} is much faster.
#'
#' @inheritParams anova_species
#'
#' @param object an object from \code{\link{dc_CA}}.
#' @param permutations a list of control values for the permutations as 
#' returned by the function \code{\link[permute]{how}}, or the number of 
#' permutations required (default 999), or a permutation matrix where each 
#' row gives the permuted indices.
#' @param rpp Logical indicating residual predictor permutation (default \code{TRUE}).
#' When \code{FALSE}, residual response permutation is used.
#' @param  by character \code{"axis"} which sets the test statistic to the 
#' first eigenvalue of the dc-CA model. Default: \code{NULL} which sets the 
#' test statistic to the inertia (sum of all double constrained eigenvalues; 
#' named  \code{constraintsTE} in the inertia element of \code{\link{dc_CA}}). 
#' This is the trait constrained inertia explained by the environmental 
#' predictors (without covariates), which is equal to the 
#' environmentally-constrained inertia explained by the traits (without trait 
#' covariates). The default is quicker computationally as it avoids
#' computation of an svd of permuted data sets.
#' 
#' @details
#' The algorithm is analogous to that of \code{\link{anova.wrda}}. The function
#' is used in \code{\link{anova.dcca}}.
#'
#' @returns
#' A list with two elements with names \code{table} and \code{eigenvalues}.
#' The \code{table} is as from \code{\link[vegan]{anova.cca}} and 
#' \code{eigenvalues} gives the dc-CA eigenvalues.
#' 
#' @references
#' ter Braak, C.J.F. & te Beest, D.E. 2022. Testing environmental effects
#' on taxonomic composition with canonical correspondence analysis:
#' alternative permutation tests are not equal.
#' Environmental and Ecological Statistics. 29 (4), 849-868.
#' \doi{10.1007/s10651-022-00545-4}
#'
#' ter Braak, C.J.F. (2022) Predictor versus response permutation
#' for significance testing in weighted regression and redundancy analysis.
#' Journal of statistical computation and simulation, 92, 2041-2059.
#' \doi{10.1080/00949655.2021.2019256}
#' 
#' @example demo/dune_test.R
#' 
#' @export
anova_sites <- function(object, 
                        permutations = 999,
                        rpp = TRUE,
                        n_axes = "all",
                        by = NULL) { 
  if (is.null(object$CWMs_orthonormal_traits)) {
    warning("Site level anova requires abundance data or ", 
            "community weighted means (CWMs).\n")
    return(list(eigenvalues = object$eigenvalues))
  }
  if (is.null(by)) by <- "omnibus"
  if (is.na(pmatch(by, c("axis", "omnibus")))) {
    stop("Set argument 'by' to 'axis' or 'NULL'.\n")
  }
  N <- nrow(object$data$dataEnv) 
  if (inherits(permutations, c("numeric", "how", "matrix"))) {
    if (is.numeric(permutations) && !is.matrix(permutations)) {
      permutations <- permute::how(nperm = permutations[1])
    } else if (is.matrix(permutations) && ncol(permutations) != N) {
      stop("Each row of permutations should have", N, "elements.\n")
    }
  } else {
    stop("Argument permutations should be integer, matrix ", 
         "or specified by permute::how().\n")
  }
  # start of new dc-ca 
  out1 <- object[c("CCAonTraits", "formulaTraits", "data", "weights", "call",
                   "Nobs", "CWMs_orthonormal_traits")]
  n <- out1$Nobs
  if (is.null(out1$CWMs_orthonormal_traits)) {
    CWMs_orthonormal_traits <- 
      scores(object$CCAonTraits, display = "species", scaling = "species", 
             choices = 1:rank_mod(object$CCAonTraits))
  } else {
    CWMs_orthonormal_traits <- out1$CWMs_orthonormal_traits * sqrt(n / (n - 1))
  }
  if (rownames(CWMs_orthonormal_traits)[1] == "col1") {
    rownames(CWMs_orthonormal_traits) <- 
      paste0("Sam", seq_len(nrow(out1$data$dataEnv)))
  }
  # step 2 Perform a weighted RDAR(M^*~E): an RDA of M^* on the 
  #        environmental variables using row weights R.
  sWn <- sqrt(object$weights$rows)
  Yw <-  CWMs_orthonormal_traits * sWn
  msqr <- msdvif(object$formulaEnv, object$data$dataEnv, object$weights$rows)
  Zw <- msqr$Zw
  Xw <- msqr$Xw
  dfpartial <- msqr$qrZ$rank
  if (rpp) {
    randperm <- randperm_eX0sqrtw
    perm_meth <- "Residualized predictor permutation\n"
  } else {
    randperm <- randperm_eY2
    perm_meth <-  "Residualized response permutation\n"
  }
  out_tes <- list()
  out_tes[[1]] <- randperm(Yw, Xw, Zw, sWn = sWn, 
                           permutations = permutations, by = by, 
                           n_axes = n_axes, return = "all")
  if (by == "axis") {
    while (out_tes[[1]]$rank > length(out_tes)) {
      Zw <- cbind(Zw, out_tes[[length(out_tes)]]$EigVector1)
      out_tes[[length(out_tes) + 1]] <- 
        randperm(Yw, Xw, Zw, sWn = sWn, permutations = permutations, 
                 by = by, return = "all")
    }
  }
  f_sites <- fanovatable(out_tes, Nobs = N, dfpartial = dfpartial, type = "row", 
                         calltext = c(object$call), perm_meth = perm_meth)
  result <- list(table = f_sites, eigenvalues = attr(f_sites, "eig"))
  return(result)
}