From d3ea03370cbaa240f370b3addd6d648b16224566 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 9 Oct 2017 16:27:30 +0300 Subject: [PATCH 01/16] add a function to amend dbrda() result with species scores Species scores are missing in dbrda, but this function can be used to add the species scores item to a dbrda() result object similarly as in other constrained ordination methods. This version is the one I suggested in https://stackoverflow.com/questions/46531969 and is still preliminary and incomplete (and undocumented and non-exported) --- R/specscores.R | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 R/specscores.R diff --git a/R/specscores.R b/R/specscores.R new file mode 100644 index 000000000..0f4dc8fac --- /dev/null +++ b/R/specscores.R @@ -0,0 +1,39 @@ +#' Add Species Scores to Ordination Results +#' +#' @param object Ordination object +#' @param comm Community data +#' +#' @rdname specscores +#' @export +`specscores` <- + function(object, comm) +{ + UseMethod("specscores") +} +#' importFrom stat +#' +#' @rdname specscores +#' @export +`specscores.dbrda` <- + function(object, comm) +{ + comm <- scale(comm, center = TRUE, scale = FALSE) + if (!is.null(object$pCCA) && object$pCCA$rank > 0) { + comm <- qr.resid(object$pCCA$QR, comm) + } + if (!is.null(object$CCA) && object$CCA$rank > 0) { + v <- crosspord(comm, object$CCA$u) + v <- decostand(v, "normalize", MARGIN = 2) + object$CCA$v <- v + comm <- qr.resid(object$CCA$QR, comm) + } + if (!is.null(object$CA) && object$CA$rank > 0) { + v <- crossprod(comm, object$CA$u) + v <- decostand(v, "normalize", MARGIN = 2) + object$CA$v <- v + } + object +} + + + From 7d6427a3f3eee20e8a4ca27f9a5061225e1aa52c Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Tue, 10 Oct 2017 14:33:57 +0300 Subject: [PATCH 02/16] species sd's needed for plot(..., correlation = TRUE) + a typo --- R/specscores.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/specscores.R b/R/specscores.R index 0f4dc8fac..b866f1879 100644 --- a/R/specscores.R +++ b/R/specscores.R @@ -18,11 +18,12 @@ function(object, comm) { comm <- scale(comm, center = TRUE, scale = FALSE) + object$colsum <- apply(comm, 2, sd) if (!is.null(object$pCCA) && object$pCCA$rank > 0) { comm <- qr.resid(object$pCCA$QR, comm) } if (!is.null(object$CCA) && object$CCA$rank > 0) { - v <- crosspord(comm, object$CCA$u) + v <- crossprod(comm, object$CCA$u) v <- decostand(v, "normalize", MARGIN = 2) object$CCA$v <- v comm <- qr.resid(object$CCA$QR, comm) From cf772cde3b4930b155d2489438d2c87acf01e069 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Tue, 10 Oct 2017 14:46:23 +0300 Subject: [PATCH 03/16] add species scores method to capscale() --- R/specscores.R | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/R/specscores.R b/R/specscores.R index b866f1879..c823852e4 100644 --- a/R/specscores.R +++ b/R/specscores.R @@ -36,5 +36,14 @@ object } +## capscale may have species scores, but is otherwise similar to dbrda + +`specscores.capscale` <- + function(object, comm) +{ + if (any(!is.na(object$colsum))) + warning("function overwrites old species scores") + specscores.dbrda(object, comm) +} From a6710e2f2e57ed5ceaf6ede93f0c7d9f544e4e17 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Tue, 10 Oct 2017 15:04:20 +0300 Subject: [PATCH 04/16] add species scores function for metaMDS Together with previous commits in this branch, this implements the functions outlined in github issue #254 for adding species scores to ordination results. Still undocumented and non-exported. --- R/specscores.R | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/R/specscores.R b/R/specscores.R index c823852e4..5541c43c5 100644 --- a/R/specscores.R +++ b/R/specscores.R @@ -46,4 +46,14 @@ specscores.dbrda(object, comm) } +## metaMDS +`specscores.metaMDS` <- + function(object, comm, expand = TRUE) +{ + if (any(!is.na(object$species))) + warning("function overwrites old species scores") + wa <- wascores(object$points, comm, expand = expand) + object$species <- wa + object +} From 019c637ced139832cf30723d4efd897ce8dcefa9 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 09:34:34 +0300 Subject: [PATCH 05/16] use replacement function to add species scores & rename to sppscores --- R/specscores.R | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/R/specscores.R b/R/specscores.R index 5541c43c5..008216254 100644 --- a/R/specscores.R +++ b/R/specscores.R @@ -5,31 +5,31 @@ #' #' @rdname specscores #' @export -`specscores` <- - function(object, comm) +`sppscores<-` <- + function(object, value) { - UseMethod("specscores") + UseMethod("sppscores<-") } #' importFrom stat #' #' @rdname specscores #' @export -`specscores.dbrda` <- - function(object, comm) +`sppscores<-.dbrda` <- + function(object, value) { - comm <- scale(comm, center = TRUE, scale = FALSE) - object$colsum <- apply(comm, 2, sd) + value <- scale(value, center = TRUE, scale = FALSE) + object$colsum <- apply(value, 2, sd) if (!is.null(object$pCCA) && object$pCCA$rank > 0) { - comm <- qr.resid(object$pCCA$QR, comm) + comm <- qr.resid(object$pCCA$QR, value) } if (!is.null(object$CCA) && object$CCA$rank > 0) { - v <- crossprod(comm, object$CCA$u) + v <- crossprod(value, object$CCA$u) v <- decostand(v, "normalize", MARGIN = 2) object$CCA$v <- v - comm <- qr.resid(object$CCA$QR, comm) + value <- qr.resid(object$CCA$QR, value) } if (!is.null(object$CA) && object$CA$rank > 0) { - v <- crossprod(comm, object$CA$u) + v <- crossprod(value, object$CA$u) v <- decostand(v, "normalize", MARGIN = 2) object$CA$v <- v } @@ -38,22 +38,22 @@ ## capscale may have species scores, but is otherwise similar to dbrda -`specscores.capscale` <- - function(object, comm) +`sppscores<-.capscale` <- + function(object, value) { if (any(!is.na(object$colsum))) warning("function overwrites old species scores") - specscores.dbrda(object, comm) + `sppscores<-.dbrda`(object, value) } ## metaMDS -`specscores.metaMDS` <- - function(object, comm, expand = TRUE) +`sppscores<-.metaMDS` <- + function(object, value) { if (any(!is.na(object$species))) warning("function overwrites old species scores") - wa <- wascores(object$points, comm, expand = expand) + wa <- wascores(object$points, value, expand = TRUE) object$species <- wa object } From 06235fd3434b86b10b181df4769938a003069519 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 12:21:45 +0300 Subject: [PATCH 06/16] renamed: R/specscores.R -> R/sppscores.R --- R/{specscores.R => sppscores.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{specscores.R => sppscores.R} (100%) diff --git a/R/specscores.R b/R/sppscores.R similarity index 100% rename from R/specscores.R rename to R/sppscores.R From 9330d2599ab6e54678fa9b0a65c04e914a077f01 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 12:24:04 +0300 Subject: [PATCH 07/16] update comments --- R/sppscores.R | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/R/sppscores.R b/R/sppscores.R index 008216254..8f62c9c84 100644 --- a/R/sppscores.R +++ b/R/sppscores.R @@ -1,19 +1,17 @@ #' Add Species Scores to Ordination Results #' #' @param object Ordination object -#' @param comm Community data +#' @param value Community data #' -#' @rdname specscores -#' @export + `sppscores<-` <- function(object, value) { UseMethod("sppscores<-") } -#' importFrom stat -#' -#' @rdname specscores -#' @export + +## dbrda + `sppscores<-.dbrda` <- function(object, value) { From 9cddaf99a1d628e271ce75de28280d63438ee03f Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 12:25:02 +0300 Subject: [PATCH 08/16] do not warn if replacement function replaces scores --- R/sppscores.R | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/R/sppscores.R b/R/sppscores.R index 8f62c9c84..cbd57d6e8 100644 --- a/R/sppscores.R +++ b/R/sppscores.R @@ -39,8 +39,6 @@ `sppscores<-.capscale` <- function(object, value) { - if (any(!is.na(object$colsum))) - warning("function overwrites old species scores") `sppscores<-.dbrda`(object, value) } @@ -49,9 +47,6 @@ `sppscores<-.metaMDS` <- function(object, value) { - if (any(!is.na(object$species))) - warning("function overwrites old species scores") - wa <- wascores(object$points, value, expand = TRUE) - object$species <- wa + object$species <- wascores(object$points, value, expand = TRUE) object } From 13ff5d061657e8ab0b856ce625ff2048f4f77b9b Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 12:52:59 +0300 Subject: [PATCH 09/16] update name of data when species scores were replaced --- R/print.metaMDS.R | 15 +++++++++------ R/sppscores.R | 4 +++- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/R/print.metaMDS.R b/R/print.metaMDS.R index 55a5df0b1..4128d74f9 100644 --- a/R/print.metaMDS.R +++ b/R/print.metaMDS.R @@ -1,5 +1,5 @@ `print.metaMDS` <- - function (x, ...) + function (x, ...) { cat("\nCall:\n") cat(deparse(x$call), "\n\n") @@ -17,29 +17,32 @@ cat(", ", c("weak", "strong")[x$ities], " ties", sep = "") cat("\n") } - if (x$converged) { - cat("Two convergent solutions found after", x$tries, + if (x$converged) { + cat("Two convergent solutions found after", x$tries, "tries\n") } else { - cat("No convergent solutions - best solution after", + cat("No convergent solutions - best solution after", x$tries, "tries\n") } z <- x$points scal <- c(if (attr(z, "centre")) "centring", if (attr(z, "pc")) "PC rotation", if (attr(z, "halfchange")) "halfchange scaling") - if (!length(scal)) + if (!length(scal)) scal <- "as is" cat("Scaling:", paste(scal, collapse = ", "), "\n") if (all(is.na(x$species))) { cat("Species: scores missing\n") } else { spattr <- attr(x$species, "shrinkage") + spdata <- attr(x$species, "data") + if (is.null(spdata)) + spdata <- x$data if (is.null(spattr)) cat("Species: non-expanded scores ") else cat("Species: expanded scores ") - cat("based on", sQuote(x$data), "\n") + cat("based on", sQuote(spdata), "\n") } cat("\n") invisible(x) diff --git a/R/sppscores.R b/R/sppscores.R index cbd57d6e8..a16e4f899 100644 --- a/R/sppscores.R +++ b/R/sppscores.R @@ -47,6 +47,8 @@ `sppscores<-.metaMDS` <- function(object, value) { - object$species <- wascores(object$points, value, expand = TRUE) + wa <- wascores(object$points, value, expand = TRUE) + attr(wa, "data") <- deparse(substitute(value)) + object$species <- wa object } From 12f6c91e6fb289730bf65b0ab547f4839990826a Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 15:45:47 +0300 Subject: [PATCH 10/16] capscale saves & print information on the data used for species scores --- R/capscale.R | 7 +++++++ R/print.cca.R | 3 +++ 2 files changed, 10 insertions(+) diff --git a/R/capscale.R b/R/capscale.R index 895a00ddc..65790d0aa 100644 --- a/R/capscale.R +++ b/R/capscale.R @@ -25,6 +25,7 @@ X <- as.dist(X) if (!inherits(X, "dist")) { comm <- X + vdata <- as.character(formula[[2]]) dfun <- match.fun(dfun) if (metaMDSdist) { commname <- as.character(formula[[2]]) @@ -35,6 +36,11 @@ } else { X <- dfun(X, distance) } + } else { # vdata name + if (missing(comm)) + vdata <- NULL + else + vdata <- deparse(substitute(comm)) } inertia <- attr(X, "method") if (is.null(inertia)) @@ -100,6 +106,7 @@ sol$CA$imaginary.u.eig <- X$negaxes } if (!is.null(comm)) { + sol$vdata <- vdata comm <- scale(comm, center = TRUE, scale = FALSE) sol$colsum <- apply(comm, 2, sd) ## take a 'subset' of the community after scale() diff --git a/R/print.cca.R b/R/print.cca.R index ef87c9a1e..685ccca5a 100644 --- a/R/print.cca.R +++ b/R/print.cca.R @@ -32,6 +32,9 @@ cs <- which(colnames(tbl) == "Rank") - 1 printCoefmat(tbl, digits = digits, na.print = "", cs.ind = seq_len(cs)) cat("Inertia is", x$inertia, "\n") + ## data used for species scores in db ordination + if (!is.null(x$vdata)) + cat("Species scores projected from", sQuote(x$vdata), "\n") if (!is.null(x$CCA$alias)) cat("Some constraints were aliased because they were collinear (redundant)\n") ## Report removed observations and species From ec32684c9b4338e9a53518493f4955da8a22ce20 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Fri, 13 Oct 2017 15:50:36 +0300 Subject: [PATCH 11/16] set the name of data used to project species scores --- R/sppscores.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/sppscores.R b/R/sppscores.R index a16e4f899..e27b3abba 100644 --- a/R/sppscores.R +++ b/R/sppscores.R @@ -15,6 +15,7 @@ `sppscores<-.dbrda` <- function(object, value) { + object$vdata <- deparse(substitute(value)) value <- scale(value, center = TRUE, scale = FALSE) object$colsum <- apply(value, 2, sd) if (!is.null(object$pCCA) && object$pCCA$rank > 0) { @@ -39,7 +40,9 @@ `sppscores<-.capscale` <- function(object, value) { - `sppscores<-.dbrda`(object, value) + object <- `sppscores<-.dbrda`(object, value) + object$vdata <- deparse(substitute(value)) + object } ## metaMDS From 5ab72a1587ea5494b9f1c5fb1ac12d590b0a4bf0 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Sat, 14 Oct 2017 14:52:13 +0300 Subject: [PATCH 12/16] update cca.object documentation: added 'vdata', no 'v' in dbrda --- man/cca.object.Rd | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/man/cca.object.Rd b/man/cca.object.Rd index b7125f88d..6b55dfdfe 100644 --- a/man/cca.object.Rd +++ b/man/cca.object.Rd @@ -195,7 +195,10 @@ \pkg{vegan} has two functions for distance-based Redundancy analysis: \code{\link{capscale}} and \code{\link{dbrda}}. Function \code{\link{capscale}} uses \code{\link{rda}} and returns its result - object, but it may add some items depending on its arguments: + object. Species scores are not directly available in distance-based + ordination but they are calculated after the analysis if community + data were available. The function can add some items depending on + its arguments: \describe{ \item{\code{metaMDSdist}}{The data set name if @@ -210,6 +213,10 @@ of sum of squares. For \eqn{n} observations this is \eqn{\sqrt{n-1}}{sqrt(n-1)} or \eqn{1} if adjustment was not done.} + \item{\code{vdata}}{Species scores (matrix \code{v}) may be + missing, but if they are present, this item gives the name of + the data that were used to calculate them.} + } Function \code{\link{dbrda}} does not use \code{\link{rda}} but @@ -225,7 +232,10 @@ eigenvalues) is given in item \code{poseig}. There is no \code{imaginary.u.eig}. The \code{ordConstrained} function returns all eigen vectors in one matrix, and the split is only made in the - \code{\link{dbrda}} function that calls \code{ordConstrained}.} } + \code{\link{dbrda}} function that calls \code{ordConstrained}.} + + \item{\code{v}}{or species scores are not calculated.} + } } From 1cf87e017e7fd87ba40b94d6bc905f4d4d39169b Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Sat, 14 Oct 2017 15:36:20 +0300 Subject: [PATCH 13/16] update documentation of capscale & dbrda This should be cherry-picked to the master even if this branch is not merged. This commit makes no reference to sppscores() development in this branch. --- man/capscale.Rd | 67 ++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/man/capscale.Rd b/man/capscale.Rd index cc266adfd..61be099a7 100644 --- a/man/capscale.Rd +++ b/man/capscale.Rd @@ -40,8 +40,8 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, they can be transformed within the formula, and they can have interactions as in a typical \code{\link{formula}}. The RHS can have a special term \code{Condition} that defines variables to be - ``partialled out'' before constraints, just like in \code{\link{rda}} - or \code{\link{cca}}. This allows the use of partial CAP.} + \dQuote{partialled out} before constraints, just like in \code{\link{rda}} + or \code{\link{cca}}. This allows the use of partial dbRDA.} \item{data}{ Data frame containing the variables on the right hand side of the model formula. } \item{distance}{The name of the dissimilarity (or distance) index if @@ -53,13 +53,15 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, \item{comm}{ Community data frame which will be used for finding species scores when the LHS of the \code{formula} was a dissimilarity matrix. This is not used if the LHS is a data - frame. If this is not supplied, the ``species scores'' are - unavailable. N.B., this is only available in \code{capscale}.} + frame. If this is not supplied, the \dQuote{species scores} are + unavailable when dissimilarities were supplied. N.B., this is + only available in \code{capscale}: \code{dbrda} does not return + species scores.} \item{add}{Add a constant to the non-diagonal dissimilarities such that all eigenvalues are non-negative in the underlying Principal Co-ordinates Analysis (see \code{\link{wcmdscale}} for - details). Choice \code{"lingoes"} (or \code{TRUE}) use the + details). \code{"lingoes"} (or \code{TRUE}) uses the recommended method of Legendre & Anderson (1999: \dQuote{method 1}) and \code{"cailliez"} uses their \dQuote{method 2}. The latter is the only one in \code{\link{cmdscale}}.} @@ -85,8 +87,8 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, expression which can contain variables in the working environment, \code{data} or species names of the community data (if given in the formula or as \code{comm} argument).} - \item{\dots}{Other parameters passed to \code{\link{rda}} or to - \code{\link{metaMDSdist}}. } + \item{\dots}{Other parameters passed to underlying functions (e.g., + \code{\link{metaMDSdist}}). } } \details{ @@ -94,7 +96,7 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, implementations of dbRDA. Function \code{capscale} is based on Legendre & Anderson (1999): the dissimilarity data are first ordinated using metric scaling, and the ordination results are - analysed with \code{\link{rda}}. Function \code{dbrda} is based on + analysed as \code{\link{rda}}. Function \code{dbrda} is based on McArdle & Anderson (2001) and directly decomposes dissimilarities. It does not use \code{\link{rda}} but a parallel implementation adapted for analysing dissimilarities and returns a @@ -111,28 +113,24 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, \code{\link{vegdist}} or distance function given in \code{dfun} with specified \code{distance}. The functions will accept distance objects from \code{\link{vegdist}}, \code{\link{dist}}, or any other - method producing similar objects. The constraining variables can be + method producing compatible objects. The constraining variables can be continuous or factors or both, they can have interaction terms, or they can be transformed in the call. Moreover, there can be a special term \code{Condition} just like in \code{\link{rda}} and - \code{\link{cca}} so that ``partial'' analysis can be performed. + \code{\link{cca}} so that \dQuote{partial} analysis can be performed. Non-Euclidean dissimilarities can produce negative eigenvalues - (Legendre & Anderson 1999, McArdle & Anderson 2001). The total - inertia and \code{\link{anova.cca}} tests for constraints will also - include the effects of imaginary axes with negative eigenvalues - following McArdle & Anderson (2001). If there are negative - eigenvalues, the printed output of \code{capscale} will add a column - with sums of positive eigenvalues and an item of sum of negative - eigenvalues, and \code{dbrda} will add a column giving the number of - real dimensions with postive eigenvalues. If negative eigenvalues - are disturbing, functions let you to distort the + (Legendre & Anderson 1999, McArdle & Anderson 2001). If there are + negative eigenvalues, the printed output of \code{capscale} will add + a column with sums of positive eigenvalues and an item of sum of + negative eigenvalues, and \code{dbrda} will add a column giving the + number of real dimensions with postive eigenvalues. If negative + eigenvalues are disturbing, functions let you to distort the dissimilarities so that only non-negative eigenvalues will be - produced with argument \code{add = TRUE} (this argument is passed - to \code{\link{cmdscale}}). Alternatively, with - \code{sqrt.dist = TRUE}, square roots of dissimilarities will be used - which may help in avoiding negative eigenvalues (Legendre & Anderson - 1999). + produced with argument \code{add = TRUE}. Alternatively, with + \code{sqrt.dist = TRUE}, square roots of dissimilarities will be + used which may help in avoiding negative eigenvalues (Legendre & + Anderson 1999). The functions can be also used to perform ordinary metric scaling a.k.a. principal coordinates analysis by using a formula with only a @@ -177,20 +175,21 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, 2003), but these developments made it similar to dbRDA. However, it discards the imaginary dimensions with negative eigenvalues and ordination and significance tests area only based on real dimensions - and positive eigenvalus. + and positive eigenvalues. The inertia is named after the dissimilarity index as defined in the - dissimilarity data, or as \code{unknown distance} if such an + dissimilarity data, or as \code{unknown distance} if such information is missing. If the largest original dissimilarity was larger than 4, \code{capscale} handles input similarly as \code{rda} - and bases its analysis on variance instead of sum of squares. Keyword - \code{mean} is added to the inertia in these cases, e.g. with - Euclidean and Manhattan distances. Inertia is based on squared index, - and keyword \code{squared} is added to the name of distance, unless - data were square root transformed (argument \code{sqrt.dist=TRUE}). If - an additive constant was used with argument \code{add}, \code{Lingoes} - or \code{Cailliez adjusted} is added to the the name of inertia, and - the value of the constant is printed.} + and bases its analysis on variance instead of sum of + squares. Keyword \code{mean} is added to the inertia in these cases, + e.g. with Euclidean and Manhattan distances. Inertia is based on + squared index, and keyword \code{squared} is added to the name of + distance, unless data were square root transformed (argument + \code{sqrt.dist=TRUE}). If an additive constant was used with + argument \code{add}, \code{Lingoes} or \code{Cailliez adjusted} is + added to the the name of inertia, and the value of the constant is + printed.} \seealso{\code{\link{rda}}, \code{\link{cca}}, \code{\link{plot.cca}}, From aadadb022afac047cedee558ecae0ac8c084dce2 Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Sun, 15 Oct 2017 18:27:58 +0300 Subject: [PATCH 14/16] export and document sppscores() --- NAMESPACE | 6 ++- man/sppscores.Rd | 108 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 113 insertions(+), 1 deletion(-) create mode 100644 man/sppscores.Rd diff --git a/NAMESPACE b/NAMESPACE index 54dc191df..c416b3aef 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,7 +27,7 @@ poolaccum, postMDS, prc, prestondistr, prestonfit, procrustes, protest, radfit, radlattice, rankindex, rarefy, rarecurve, rareslope, raupcrick, rda, renyiaccum, renyi, rrarefy, scores, showvarparts, simper, spandepth, spantree, specaccum, specslope, -specnumber, specpool2vect, specpool, spenvcor, +specnumber, specpool2vect, specpool, spenvcor, "sppscores<-", stepacross, stressplot, swan, tabasco, taxa2dist, taxondive, tolerance, treedist, treedive, treeheight, tsallisaccum, tsallis, varpart, vectorfit, vegandocs, vegdist, avgdist, vegemite, veiledspec, wascores, @@ -467,6 +467,10 @@ S3method(simulate, nullmodel) # specslope: vegan S3method(specslope, specaccum) S3method(specslope, fitspecaccum) +## sppscores<-: vegan +S3method("sppscores<-", capscale) +S3method("sppscores<-", dbrda) +S3method("sppscores<-", metaMDS) # SSD: stats S3method(SSD, cca) # str: utils diff --git a/man/sppscores.Rd b/man/sppscores.Rd new file mode 100644 index 000000000..95d6f77e7 --- /dev/null +++ b/man/sppscores.Rd @@ -0,0 +1,108 @@ +\name{sppscores} +\alias{sppscores<-} +\alias{sppscores<-.dbrda} +\alias{sppscores<-.capscale} +\alias{sppscores<-.metaMDS} + +\title{ +Add or Replace Species Scores in Distance-Based Ordination +} + +\description{ + + Distance-based ordination (\code{\link{dbrda}}, + \code{\link{capscale}}, \code{\link{metaMDS}}) have no information + on species, but some methods may add species scores if community + data were available. However, the species scores may be missing (and + they always are in \code{\link{dbrda}}), or they may not have a + close relation to used dissimilarity index. This function will add + the species scores or replace the existing species scores in + distance-based methods. + +} + +\usage{ +sppscores(object) <- value +} + +\arguments{ + \item{object}{Ordination result.} + \item{value}{Community data to find the species scores.} +} + +\details{ + + Distances have no information on species (columns, variables), and + hence distance-based ordination has no information on species + scores. However, the species scores can be added as supplementary + information after the analysis to help the interpretation of + results. Some ordination methods (\code{\link{capscale}}, + \code{\link{metaMDS}}) can supplement the species scores during the + analysis if community data was available in the analysis. + + In \code{\link{capscale}} the species scores are found by projecting + the community data to site ordination (linear combination scores), + and the scores are accurate if the analysis used Euclidean + distances. If the dissimilarity index can be expressed as Euclidean + distances of transformed data (for instance, Chord and Hellinger + Distances), the species scores based on transformed data will be + accurate, but the function still finds the dissimilarities with + untransformed data. Usually community dissimilarities differ in two + significant ways from Euclidean distances: They are bound to maximum + 1, and they use absolute differences instead of squared + differences. In such cases, it may be better to use species scores + that are transformed so that their Euclidean distances have a good + linear relation to used dissimilarities. It is often useful to + standardize data so that each row has unit total, and perform + squareroot transformation to damp down the effect of squared + differences (see Examples). + + Function \code{\link{dbrda}} never finds the species scores, but it + is mathematically similar to \code{\link{capscale}}, and similar + rules should be followed when supplementing the species scores. + + Function \code{\link{metaMDS}} uses weighted averages + (\code{\link{wascores}}) to find the species scores. These have a + better relationship with most dissimilarities than the projection + scores used in metric ordination, but similar transformation of the + community data should be used both in dissimilarities and in species + scores. + +} + +\value{ + + Replacement function adds the species scores or replaces the old + scores in the ordination object. + +} + + +\author{ + Jari Oksanen +} + +\seealso{ + + Function \code{\link{envfit}} finds similar scores, but based on + correlations. The species scores for non-metric ordination use + \code{\link{wascores}} which can also used directly on any + ordination result. + +} + +\examples{ +data(BCI, BCI.env) +mod <- dbrda(vegdist(BCI) ~ Habitat, BCI.env) +## add species scores +sppscores(mod) <- BCI +## Euclidean distances of BCI differ from used dissimilarity +plot(vegdist(BCI), dist(BCI)) +## more linear relationship +plot(vegdist(BCI), dist(sqrt(decostand(BCI, "total")))) +## better species scores +sppscores(mod) <- sqrt(decostand(BCI, "total")) +} + +\keyword{ multivariate } + From fa1ff35ba1ce68ca3da020575d74e7a941e15c1c Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 16 Oct 2017 09:49:54 +0300 Subject: [PATCH 15/16] token implementation of sppscores() --- NAMESPACE | 2 +- R/sppscores.R | 10 ++++++++++ man/sppscores.Rd | 1 + 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/NAMESPACE b/NAMESPACE index c416b3aef..b464e47df 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,7 +27,7 @@ poolaccum, postMDS, prc, prestondistr, prestonfit, procrustes, protest, radfit, radlattice, rankindex, rarefy, rarecurve, rareslope, raupcrick, rda, renyiaccum, renyi, rrarefy, scores, showvarparts, simper, spandepth, spantree, specaccum, specslope, -specnumber, specpool2vect, specpool, spenvcor, "sppscores<-", +specnumber, specpool2vect, specpool, spenvcor, "sppscores<-", sppscores, stepacross, stressplot, swan, tabasco, taxa2dist, taxondive, tolerance, treedist, treedive, treeheight, tsallisaccum, tsallis, varpart, vectorfit, vegandocs, vegdist, avgdist, vegemite, veiledspec, wascores, diff --git a/R/sppscores.R b/R/sppscores.R index e27b3abba..e80fff7cc 100644 --- a/R/sppscores.R +++ b/R/sppscores.R @@ -55,3 +55,13 @@ object$species <- wa object } + +## the main purpose of accessor function is to provide nicer command +## autocompletion and cross-references in help, and of course, to tell +## that it is not implemented (and may never be) + +`sppscores` <- + function(object) +{ + .NotYetImplemented() +} diff --git a/man/sppscores.Rd b/man/sppscores.Rd index 95d6f77e7..79cb3ac87 100644 --- a/man/sppscores.Rd +++ b/man/sppscores.Rd @@ -1,4 +1,5 @@ \name{sppscores} +\alias{sppscores} \alias{sppscores<-} \alias{sppscores<-.dbrda} \alias{sppscores<-.capscale} From d11f57113f9143d2708fa824268531935f5fdbed Mon Sep 17 00:00:00 2001 From: Jari Oksanen Date: Mon, 16 Oct 2017 10:05:16 +0300 Subject: [PATCH 16/16] dbrda/capscale documentation tells about sppscores() --- man/capscale.Rd | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/man/capscale.Rd b/man/capscale.Rd index 61be099a7..399fa88fa 100644 --- a/man/capscale.Rd +++ b/man/capscale.Rd @@ -56,7 +56,8 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, frame. If this is not supplied, the \dQuote{species scores} are unavailable when dissimilarities were supplied. N.B., this is only available in \code{capscale}: \code{dbrda} does not return - species scores.} + species scores. Function \code{\link{sppscores}} can be used to add + species scores if they are missing.} \item{add}{Add a constant to the non-diagonal dissimilarities such that all eigenvalues are non-negative in the underlying Principal @@ -119,6 +120,10 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, special term \code{Condition} just like in \code{\link{rda}} and \code{\link{cca}} so that \dQuote{partial} analysis can be performed. + Function \code{dbrda} does not return species scores, and they can + also be missing in \code{capscale}, but they can be added after the + analysis using function \code{\link{sppscores}}. + Non-Euclidean dissimilarities can produce negative eigenvalues (Legendre & Anderson 1999, McArdle & Anderson 2001). If there are negative eigenvalues, the printed output of \code{capscale} will add @@ -194,7 +199,9 @@ dbrda(formula, data, distance = "euclidean", sqrt.dist = FALSE, \seealso{\code{\link{rda}}, \code{\link{cca}}, \code{\link{plot.cca}}, \code{\link{anova.cca}}, \code{\link{vegdist}}, - \code{\link{dist}}, \code{\link{cmdscale}}, \code{\link{wcmdscale}}. + \code{\link{dist}}, \code{\link{cmdscale}}, \code{\link{wcmdscale}} + for underlying and related functions. Function \code{\link{sppscores}} + can add species scores or replace existing species scores. The function returns similar result object as \code{\link{rda}} (see \code{\link{cca.object}}). This section for \code{\link{rda}} gives a