[R-sig-eco] CCA outputs in Vegan

[Gavin Simpson]

On Fri, 2011-02-25 at 15:55 -0400, Dominique wrote:
> Dear all,
> 
>  
> 
> I am desperately searching for code lines able to modify the standard
> outputs of CCA in Vegan. I need to have an access to the three
> parameters U, C and T, usually employed for unimodal species response
> curves:
> 
> U : centre of species distributions (species scores) is presently
> given but
> 
> C : species maximal theoretical abundance, named  N2  in 'old' Canoco
> outputs and

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.

> T : "species tolerance" on each axis = (the range of species occurrences is 4T).

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

>  
> 
> ..are lacking
> 
>  
> 
> It would be helpful if anyone knew the syntax required
> 
>  
> 
> Best regards
> 
> Dominique
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
-------------- next part --------------
#' 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 choices numeric; the ordination scaling to use.
#' @param useN2 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)
}