[Vegan-commits] r2370 - in pkg/vegan: R inst

[noreply at r-forge.r-project.org]

Author: jarioksa
Date: 2013-01-15 09:13:02 +0100 (Tue, 15 Jan 2013)
New Revision: 2370

Modified:
   pkg/vegan/R/stressplot.wcmdscale.R
   pkg/vegan/inst/ChangeLog
Log:
stressplot.capscale also uses dist(u.eig %*% t(v)) instead of dist(u.eig)

Modified: pkg/vegan/R/stressplot.wcmdscale.R
===================================================================
--- pkg/vegan/R/stressplot.wcmdscale.R	2013-01-14 18:49:04 UTC (rev 2369)
+++ pkg/vegan/R/stressplot.wcmdscale.R	2013-01-15 08:13:02 UTC (rev 2370)
@@ -110,18 +110,35 @@
     ## Scores to reconstruct data
     u <- cbind(object$CCA$u, object$CA$u)
     ev <- c(object$CCA$eig, object$CA$eig)
-    u <- u %*% diag(sqrt(ev))
-    if (!is.null(object$pCCA))
+    if (object$adjust == 1)
+        const <- sqrt(NROW(u) - 1)
+    else
+        const <- 1
+    u <- u %*% diag(sqrt(ev) * const)
+    ## Constrained ordination needs also scores 'v' to reconstruct
+    ## 'data', but these are not returned by capscale() which replaces
+    ## original 'v' with weighted sums of 'comm' data.
+    if (!is.null(object$CCA)) 
+        v <- svd(object$CCA$Xbar - object$CA$Xbar, nu = 0, nv = object$CCA$qrank)$v
+    else
+        v <- NULL
+    if (!is.null(object$CA))
+        v <- cbind(v, svd(object$CA$Xbar, nu = 0, nv = object$CA$rank)$v)
+    ## Reconstruct Xbar and Xbark
+    Xbar <- u %*% t(v)
+    Xbark <- u[,seq_len(k), drop = FALSE] %*% t(v[,seq_len(k), drop = FALSE])
+    if (!is.null(object$pCCA)) {
         pFit <- object$pCCA$Fit/object$adjust
-    else
-        pFit <- NULL
+        Xbar <- Xbar + pFit
+        Xbark <- Xbark + pFit
+    }
     ## Distances
-    dis <- dist(cbind(u, pFit))
+    dis <- dist(Xbar)
+    odis <- dist(Xbark)
     if (!is.null(object$CA$imaginary.u.eig))
         dis <- sqrt(dis^2 - dist(object$CA$imaginary.u.eig)^2)
     if (!is.null(object$ac))
         dis <- dis - object$ac
-    odis <- dist(cbind(u[,seq_len(k), drop=FALSE], pFit))
     ## plot like above
         ## Plot
     if (missing(pch))

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2013-01-14 18:49:04 UTC (rev 2369)
+++ pkg/vegan/inst/ChangeLog	2013-01-15 08:13:02 UTC (rev 2370)
@@ -22,8 +22,7 @@
 
 	* stressplot: added stressplot() methods for wcmdscale(),
 	capscale(), cca(), rda(), prcomp() and princomp() results.  These
-	also work with constrained ordination, but not with the partial
-	ordination (at least not yet).  These methods display the
+	also work with constrained ordination. 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
@@ -33,6 +32,17 @@
 	models (p-dbRDA, pRDA, pCCA) add the partial component both to the
 	original dissimilarities and the fit.
 
+	The row scores (u) alone will not correctly estimate original
+	dissimilarities in constrained (or partial) ordination. In
+	unconstrained ordination we can get the distances as dist(u %*%
+	diag(sqrt(eig))), but in constrained ordination this will not give
+	the observed dissimilarities with all axes. Currently we get the
+	ordination distances from a (low-rank) approximation of the data
+	as dist(u %*% diag(sqrt(eig)) %*% t(v)). However, it is not sure
+	that this the right thing to do, but perhaps we should acknowledge
+	the fact row ordination with constraints does not approximate
+	distances. So this may change.
+
 	* wcmdscale: added method functions print(), plot() and
 	scores(). Now class "wcmdscale" results also retun the function
 	call and dimensions have names.