Issue #254

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<-", 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

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()

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

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)

R/sppscores.R

@@ -0,0 +1,67 @@
+#' Add Species Scores to Ordination Results
+#'
+#' @param object Ordination object
+#' @param value Community data
+#'
+
+`sppscores<-` <-
+    function(object, value)
+{
+    UseMethod("sppscores<-")
+}
+
+## dbrda
+
+`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) {
+        comm <- qr.resid(object$pCCA$QR, value)
+    }
+    if (!is.null(object$CCA) && object$CCA$rank > 0) {
+        v <- crossprod(value, object$CCA$u)
+        v <- decostand(v, "normalize", MARGIN = 2)
+        object$CCA$v <- v
+        value <- qr.resid(object$CCA$QR, value)
+    }
+    if (!is.null(object$CA) && object$CA$rank > 0) {
+        v <- crossprod(value, object$CA$u)
+        v <- decostand(v, "normalize", MARGIN = 2)
+        object$CA$v <- v
+    }
+    object
+}
+
+## capscale may have species scores, but is otherwise similar to dbrda
+
+`sppscores<-.capscale` <-
+    function(object, value)
+{
+    object <- `sppscores<-.dbrda`(object, value)
+    object$vdata <- deparse(substitute(value))
+    object
+}
+
+## metaMDS
+
+`sppscores<-.metaMDS` <-
+    function(object, value)
+{
+    wa <- wascores(object$points, value, expand = TRUE)
+    attr(wa, "data") <- deparse(substitute(value))
+    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()
+}

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,16 @@ 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. 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
      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,16 +88,16 @@ 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{
 
   Functions \code{capscale} and \code{dbrda} provide two alternative
   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 +114,28 @@ 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.
+
+  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). 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,25 +180,28 @@ 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}},
   \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

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.}
+   }
 
   }