Skip to content

CCA outputs in Vegan

2 messages · Dominique, Gavin Simpson

Original
2 days later
Gavin Simpson
#
On Fri, 2011-02-25 at 15:55 -0400, Dominique wrote:
N2 can be computed using `renyi()` in vegan. The attached function can
compute it for you, and return it via argument useN2 = TRUE, however
note that it will alter the tolerance / heterogeneity values, and Canoco
doesn't apply the N2 bias correction in the sol file., but Canodraw will
use it.
I was asked about the latter some time ago and made some progress before
being distracted by the day job. Your email has spurred me on to tackle
this again.

Note; I haven't tested this code against Canoco and it does not
reproduce the example shown in the Canoco manual as I'm not a hundred
per cent sure what analysis that pertains to, and it appears to have
used an augmented version of the Dune data set. I also don't have access
to Canoco easily via my laptop (running Linux) so haven't been able to
run any checks.

As fas as I can tell, the function does implement the two relevant
equations but it has been optimised to use R-isms to speed the whole
thing up. The optimisations replicated the direct translation of the two
relevant equations.

I'd appreciate some eyes on this if anyone wants to use the dune data as
provided by Vegan and run it through the same analysis in Canoco to
check the results.

Anyway, here is the function (and attached):

#' Species tolerances and sample heterogeneities
#'
#' Function to compute species tolerances and site heterogeneity measures
#' from unimodal ordinations (CCA & CA). Implements Eq 6.47 and 6.48 from
#' the Canoco 4.5 Reference Manual (pages 178-179).
#'
#' Licence: GPL Version 2 http://www.gnu.org/licenses/gpl-2.0.html
#'
#' @param x object of class \code{"cca"}.
#' @param Y original species data matrix.
#' @param choices numeric; which ordination axes to compute tolerances
#'   and heterogeneities for. Defaults to axes 1 and 2.
#' @param which; character, one of \code{"species"} or \code{"sites"},
#'   indicating whether species tolerances or sample heterogeneities
#'   respectively are computed.
#' @param numeric; the ordination scaling to use.
#' @param logical; should the bias in the tolerances/heterogeneities
#'   be reduced via scaling by Hill's N2?
#' @return matrix of tolerances/heterogeneities with some additional
#'   attributes.
#' @author Gavin Simpson \email{gavin.simpson AT ucl.ac.uk}
#' @examples
#' data(dune)
#' data(dune.env)
#' mod <- cca(dune ~ ., data = dune.env)
#' tolerance.cca(mod, dune)
tolerance.cca <- function(x, Y, choices = 1:2,
                          which = c("species","sites"),
                          scaling = 2, useN2 = FALSE) {
    if(inherits(x, "rda"))
        stop("Tolerances only available for unimodal ordinations.")
    if(missing(which))
        which <- "species"
    which <- match.arg(which)
    siteScrTypes <- if(is.null(x$CCA)){ "sites" } else {"lc"}
    scrs <- scores(x, display = c(siteScrTypes,"species"),
                   choices = choices, scaling = scaling)
    ## compute N2 if useN2 == TRUE & only if
    doN2 <- isTRUE(useN2) && ((which == "species" && abs(scaling) == 2) ||
                              (which == "sites" && abs(scaling) == 1))

    ## this gives the x_i - u_k on axis j
    ## outer(scrs$sites, scrs$species, "-")[,2,,j]
    siteScrs <- which(names(scrs) %in% c("sites","constraints"))
    xiuk <- outer(scrs[[siteScrs]], scrs$species, "-")
    if(isTRUE(all.equal(which, "sites"))) {
        ## need to permute the array as rowSums has different idea of what rows
        ## are that doesn't correspond to colSums. So flip dimensions 1 and 2
        ## with aperm and use colSums.
        res <- sqrt(sweep(colSums(aperm(sweep(xiuk[ , 2, , choices]^2, c(1:2),
                                              data.matrix(Y), "*"),
                                        c(2,1,3))),
                          1, rowSums(Y), "/"))
        if(doN2) {
            tot <- rowSums(Y)
            y <- sweep(Y, 1, tot, "/")^2
            N2 <- 1 / rowSums(y, na.rm = TRUE) ## 1/H
            res <- sweep(res, 1, sqrt(1 - (1/N2)), "/")
        }
    } else {
        res <- sqrt(sweep(colSums(sweep(xiuk[ , 2, , choices]^2, c(1:2),
                                        data.matrix(Y), "*")),
                          1, colSums(Y), "/"))
        if(doN2) {
            tot <- colSums(Y)
            y <- sweep(Y, 2, tot, "/")^2
            N2 <- 1 / colSums(y, na.rm = TRUE) ## 1/H
            res <- sweep(res, 1, sqrt(1 - (1/N2)), "/")
        }
    }
    class(res) <- c("tolerance.cca","tolerance","matrix")
    attr(res, "which") <- which
    attr(res, "scaling") <- scaling
    attr(res, "N2") <- NULL
    if(doN2)
        attr(res, "N2") <- N2
    attr(res, "model") <- deparse(substitute(mod))
    return(res)
}

HTH

G