[Vegan-commits] r2353 - in pkg/vegan: . R inst man tests/Examples

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Jan 9 11:21:14 CET 2013


Author: jarioksa
Date: 2013-01-09 11:21:14 +0100 (Wed, 09 Jan 2013)
New Revision: 2353

Added:
   pkg/vegan/R/stressplot.wcmdscale.R
   pkg/vegan/man/stressplot.wcmdscale.Rd
Modified:
   pkg/vegan/NAMESPACE
   pkg/vegan/inst/ChangeLog
   pkg/vegan/man/goodness.metaMDS.Rd
   pkg/vegan/tests/Examples/vegan-Ex.Rout.save
Log:
add stressplot methods for cca, rda, capscale & wcmdscale

Modified: pkg/vegan/NAMESPACE
===================================================================
--- pkg/vegan/NAMESPACE	2013-01-09 08:50:25 UTC (rev 2352)
+++ pkg/vegan/NAMESPACE	2013-01-09 10:21:14 UTC (rev 2353)
@@ -401,6 +401,10 @@
 # stressplot: vegan
 S3method(stressplot, default)
 S3method(stressplot, monoMDS)
+S3method(stressplot, wcmdscale)
+S3method(stressplot, capscale)
+S3method(stressplot, cca)
+S3method(stressplot, rda)
 # summary: base
 S3method(summary, anosim)
 S3method(summary, bioenv)

Added: pkg/vegan/R/stressplot.wcmdscale.R
===================================================================
--- pkg/vegan/R/stressplot.wcmdscale.R	                        (rev 0)
+++ pkg/vegan/R/stressplot.wcmdscale.R	2013-01-09 10:21:14 UTC (rev 2353)
@@ -0,0 +1,112 @@
+### stressplot() methods for eigenvector ordinations wcmdscale, rda,
+### cca, capscale
+
+`stressplot.wcmdscale` <-
+    function(object, k = 2, pch,  p.col = "blue", l.col = "red", lwd = 2, ...)
+{
+    ## NB: Ignores weights
+    
+    ## Get the ordination distances in k dimensions
+    if (k > NCOL(object$points))
+        stop("'k' cannot exceed the number of real dimensions")
+    odis <- dist(object$points[,1:k, drop = FALSE])
+    ## Reconstitute the original observed distances
+    dis <- dist(object$points)
+    if (!is.null(object$negaxes))
+        dis <- sqrt(dis^2 - dist(object$negaxes)^2)
+    ## Plot
+    if (missing(pch))
+        if (length(dis) > 5000)
+            pch <- "."
+        else
+            pch <- 1
+    plot(dis, odis, pch = pch, col = p.col, xlab = "Observed Dissimilarity",
+         ylab = "Ordination Distance", ...)
+    abline(0, 1, col = l.col, lwd = lwd, ...)
+    invisible(odis)
+}
+
+`stressplot.rda` <-
+    function(object, k = 2, pch, p.col = "blue", l.col = "red", lwd = 2, ...)
+{
+    ## Not yet done for pRDA
+    if (!is.null(object$pCCA))
+        stop("not implemented yet for partial RDA")
+    ## Normalized scores to reconstruct data
+    u <- cbind(object$CCA$u, object$CA$u)
+    ev <- c(object$CCA$eig, object$CA$eig)
+    ## normalizing constant
+    nr <- NROW(u)
+    const <- sqrt(ev * (nr-1))
+    ## Distances
+    dis <- dist(u %*% diag(const))
+    odis <- dist(u[,1:k, drop=FALSE] %*% diag(const[1:k], nrow = k))
+    ## plot like above
+        ## Plot
+    if (missing(pch))
+        if (length(dis) > 5000)
+            pch <- "."
+        else
+            pch <- 1
+    plot(dis, odis, pch = pch, col = p.col, xlab = "Observed Dissimilarity",
+         ylab = "Ordination Distance", ...)
+    abline(0, 1, col = l.col, lwd = lwd, ...)
+    invisible(odis)
+}
+
+`stressplot.cca` <-
+    function(object, k = 2, pch, p.col = "blue", l.col = "red", lwd = 2, ...)
+{
+    ## NB: Ignores row weights!
+    
+    ## Not yet done for pCCA
+    if (!is.null(object$pCCA))
+        stop("not implemented yet for partial CCA")
+    ## Normalized scores to reconstruct data
+    u <- cbind(object$CCA$u, object$CA$u)
+    sev <- sqrt(c(object$CCA$eig, object$CA$eig))
+    ## Distances
+    dis <- dist(u %*% diag(sev))
+    odis <- dist(u[,1:k, drop=FALSE] %*% diag(sev[1:k], nrow = k))
+    ##odis <- dist(sweep(Xbar, 2, sqrt(object$colsum), "*"))
+    ## plot like above
+        ## Plot
+    if (missing(pch))
+        if (length(dis) > 5000)
+            pch <- "."
+        else
+            pch <- 1
+    plot(dis, odis, pch = pch, col = p.col, xlab = "Observed Dissimilarity",
+         ylab = "Ordination Distance", ...)
+    abline(0, 1, col = l.col, lwd = lwd, ...)
+    invisible(odis)
+}
+
+`stressplot.capscale` <-
+    function(object, k = 2, pch, p.col = "blue", l.col = "red", lwd = 2, ...)
+{
+    ## Not yet done for pRDA
+    if (!is.null(object$pCCA))
+        stop("not implemented yet for partial analysis")
+    ## Normalized scores to reconstruct data
+    u <- cbind(object$CCA$u, object$CA$u)
+    ev <- c(object$CCA$eig, object$CA$eig)
+    ## normalizing constant
+    const <- sqrt(ev)
+    ## Distances
+    dis <- dist(u %*% diag(const))
+    if (!is.null(object$CA$imaginary.u.eig))
+        dis <- sqrt(dis^2 - dist(object$CA$imaginary.u.eig)^2)
+    odis <- dist(u[,1:k, drop=FALSE] %*% diag(const[1:k], nrow = k))
+    ## plot like above
+        ## Plot
+    if (missing(pch))
+        if (length(dis) > 5000)
+            pch <- "."
+        else
+            pch <- 1
+    plot(dis, odis, pch = pch, col = p.col, xlab = "Observed Dissimilarity",
+         ylab = "Ordination Distance", ...)
+    abline(0, 1, col = l.col, lwd = lwd, ...)
+    invisible(odis)
+}

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2013-01-09 08:50:25 UTC (rev 2352)
+++ pkg/vegan/inst/ChangeLog	2013-01-09 10:21:14 UTC (rev 2353)
@@ -10,6 +10,17 @@
 	with decostand() function. This change means that scaling of
 	species scores will change, and graphs can look different than
 	previously. All analyses should be redone. 
+
+	* stressplot: added stressplot() methods for wcmdscale(),
+	capscale(), cca() and rda() results.  These also work with
+	constrained ordination, but not with the partial ordination (at
+	least not yet).  These methods display the ordination distances in
+	given number of dimensions (defaults 'k = 2') against original
+	observed distances.  These original distances are found from the
+	full space solution, and in capscale() and wcmdscale() they are
+	correct for the imaginary axes. The weights are not used in
+	wcmdscale() and cca(): I must first figure out if they should be
+	used or ignored.
 	
 Version 2.1-22 (closed January 8, 2013)
 

Modified: pkg/vegan/man/goodness.metaMDS.Rd
===================================================================
--- pkg/vegan/man/goodness.metaMDS.Rd	2013-01-09 08:50:25 UTC (rev 2352)
+++ pkg/vegan/man/goodness.metaMDS.Rd	2013-01-09 10:21:14 UTC (rev 2353)
@@ -66,13 +66,17 @@
 } 
 
 \value{ Function \code{goodness} returns a vector of values. Function
-  \code{stressplot} returns invisibly an object with itmes for
+  \code{stressplot} returns invisibly an object with items for
   original dissimilarities, ordination distances and fitted values.  }
 
 \author{Jari Oksanen. }
 
 \seealso{\code{\link{metaMDS}},  \code{\link{monoMDS}}, 
-  \code{\link[MASS]{isoMDS}}, \code{\link[MASS]{Shepard}}. }
+  \code{\link[MASS]{isoMDS}}, \code{\link[MASS]{Shepard}}. Similar
+  diagrams for eigenvector ordinations can be drawn with
+  \code{\link{stressplot.wcmdscale}}, \code{\link{stressplot.cca}},
+  \code{\link{stressplot.rda}} and \code{\link{stressplot.capscale}}.
+}
 
 \examples{
 data(varespec)

Added: pkg/vegan/man/stressplot.wcmdscale.Rd
===================================================================
--- pkg/vegan/man/stressplot.wcmdscale.Rd	                        (rev 0)
+++ pkg/vegan/man/stressplot.wcmdscale.Rd	2013-01-09 10:21:14 UTC (rev 2353)
@@ -0,0 +1,81 @@
+\name{stressplot.wcmdscale}
+\alias{stressplot.wcmdscale}
+\alias{stressplot.cca}
+\alias{stressplot.rda}
+\alias{stressplot.capscale}
+
+\title{
+  Display Ordination Distances Against Observed Distances in Eigenvector Ordinations
+}
+
+\description{
+  Functions plot ordination distances in given number of dimensions
+  against observed distances or distances in full space in eigenvector
+  methods. The display is similar as the Shepard diagram
+  (\code{\link{stressplot}} for non-metric multidimensional scaling
+  with \code{\link{metaMDS}} or \code{\link{monoMDS}}), but shows the
+  linear relationship of the eigenvector ordinations. The
+  \code{stressplot} methods are available for \code{\link{wcmdscale}},
+  \code{\link{rda}}, \code{\link{cca}} and \code{\link{capscale}}. 
+}
+
+\usage{
+\method{stressplot}{wcmdscale}(object, k = 2, pch, p.col = "blue", l.col = "red",
+    lwd = 2, ...)
+}
+
+\arguments{
+  \item{object}{
+    Result object from eigenvector ordination (\code{\link{wcmdscale}},
+    \code{\link{rda}}, \code{\link{cca}}, \code{\link{capscale}})
+}
+  \item{k}{
+    Number of dimensions for which the ordination distances are displayed.
+}
+  \item{pch, p.col, l.col, lwd}{
+    Plotting character, point colour and line colour like in
+    default \code{\link{stressplot}}
+}
+  \item{\dots}{
+    Other parameters to fuctions, e.g. graphical parameters.
+}
+}
+
+\details{ The functions offer a similar display for eigenvector
+  ordinations as the standard Shepard diagram (\code{\link{stressplot}})
+  in non-metric multidimensional scaling. The ordination distances in
+  given number of dimensions are plotted against observed
+  distances. With metric distances, the ordination distances in full
+  space (with all ordination axes) are equal to observed distances, and
+  the fit line shows this equality. In general, the fit line does not go
+  through the points, but the points for observed distances approach the
+  fit line from below. However, with non-metric distances (in
+  \code{\link{wcmdscale}} or \code{\link{capscale}}) with negative
+  eigenvalues the ordination distances can exceed the observed distances
+  in real dimensions; the imaginary dimensions with negative eigenvalues
+  will correct these excess distances.
+}
+
+\value{
+  Functions draw a graph and return invisibly the ordination distances.
+}
+
+\author{
+  Jari Oksanen.
+}
+
+\seealso{
+  \code{\link{stressplot}} and \code{\link{stressplot.monoMDS}} for
+  standard Shepard diagrams.
+}
+
+\examples{
+data(dune, dune.env)
+mod <- rda(dune)
+stressplot(mod)
+mod <- rda(dune ~ Management, dune.env)
+stressplot(mod, k=3)
+}
+
+\keyword{ multivariate }
+

Modified: pkg/vegan/tests/Examples/vegan-Ex.Rout.save
===================================================================
--- pkg/vegan/tests/Examples/vegan-Ex.Rout.save	2013-01-09 08:50:25 UTC (rev 2352)
+++ pkg/vegan/tests/Examples/vegan-Ex.Rout.save	2013-01-09 10:21:14 UTC (rev 2353)
@@ -161,7 +161,7 @@
 
 Formula:
 y ~ poly(x1, 1) + poly(x2, 1)
-<environment: 0x12bef70>
+<environment: 0x3160ba8>
 Total model degrees of freedom 3 
 
 GCV score: 0.04278782
@@ -4910,7 +4910,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x8b27e80>
+<environment: 0xa20afb0>
 
 Estimated degrees of freedom:
 6.45  total = 7.45 
@@ -4926,7 +4926,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x905b218>
+<environment: 0xa251f70>
 
 Estimated degrees of freedom:
 6.12  total = 7.12 
@@ -5086,7 +5086,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x8434c00>
+<environment: 0xa336010>
 
 Estimated degrees of freedom:
 8.93  total = 9.93 
@@ -5099,7 +5099,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x98d3c30>
+<environment: 0xab88d18>
 
 Estimated degrees of freedom:
 7.75  total = 8.75 
@@ -5112,7 +5112,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x9b013d0>
+<environment: 0xab15780>
 
 Estimated degrees of freedom:
 8.9  total = 9.9 
@@ -6870,6 +6870,29 @@
 > 
 > 
 > cleanEx()
+> nameEx("stressplot.wcmdscale")
+> ### * stressplot.wcmdscale
+> 
+> flush(stderr()); flush(stdout())
+> 
+> ### Name: stressplot.wcmdscale
+> ### Title: Display Ordination Distances Against Observed Distances in
+> ###   Eigenvector Ordinations
+> ### Aliases: stressplot.wcmdscale stressplot.cca stressplot.rda
+> ###   stressplot.capscale
+> ### Keywords: multivariate
+> 
+> ### ** Examples
+> 
+> data(dune, dune.env)
+> mod <- rda(dune)
+> stressplot(mod)
+> mod <- rda(dune ~ Management, dune.env)
+> stressplot(mod, k=3)
+> 
+> 
+> 
+> cleanEx()
 > nameEx("taxondive")
 > ### * taxondive
 > 
@@ -7666,7 +7689,7 @@
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0xa84c790>
+<environment: 0x950e978>
 
 Estimated degrees of freedom:
 2  total = 3 
@@ -8146,7 +8169,7 @@
 > ### * <FOOTER>
 > ###
 > cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed:  26.141 0.1 26.342 0 0 
+Time elapsed:  26.121 0.156 26.372 0 0 
 > grDevices::dev.off()
 null device 
           1 



More information about the Vegan-commits mailing list