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

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

Author: jarioksa
Date: 2013-01-25 13:51:10 +0100 (Fri, 25 Jan 2013)
New Revision: 2388

Modified:
   pkg/vegan/DESCRIPTION
   pkg/vegan/R/protest.R
   pkg/vegan/inst/ChangeLog
Log:
huge speed-up in protest

Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION	2013-01-25 12:07:42 UTC (rev 2387)
+++ pkg/vegan/DESCRIPTION	2013-01-25 12:51:10 UTC (rev 2388)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 2.1-23
-Date: January 8, 2013
+Version: 2.1-24
+Date: January 25, 2013
 Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, 
    Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, 
    M. Henry H. Stevens, Helene Wagner  

Modified: pkg/vegan/R/protest.R
===================================================================
--- pkg/vegan/R/protest.R	2013-01-25 12:07:42 UTC (rev 2387)
+++ pkg/vegan/R/protest.R	2013-01-25 12:51:10 UTC (rev 2388)
@@ -3,9 +3,27 @@
 {
     X <- scores(X, display = scores, ...)
     Y <- scores(Y, display = scores, ...)
-    sol <- procrustes(X, Y, symmetric = TRUE)
+    ## Centre and normalize X & Y here so that the permutations will
+    ## be faster
+    X <- scale(X, scale = FALSE)
+    Y <- scale(Y, scale = FALSE)
+    X <- X/sqrt(sum(X^2))
+    Y <- Y/sqrt(sum(Y^2))
+    ## Transformed X and Y will yield symmetric procrustes() and we
+    ## need not specify that in the call (but we set it symmetric
+    ## after the call).
+    sol <- procrustes(X, Y, symmetric = FALSE)
+    sol$symmetric <- TRUE
     sol$t0 <- sqrt(1 - sol$ss)
     N <- nrow(X)
+
+    ## Permutations: We only need the goodness of fit statistic from
+    ## Procrustes analysis, and therefore we only have the necessary
+    ## function here. This avoids a lot of overhead of calling
+    ## procrustes() for each permutation. The following gives the
+    ## Procrustes r directly.
+    procr <- function(X, Y) sum(svd(crossprod(X, Y), nv=0, nu=0)$d)
+    
     if (length(permutations) == 1) {
         if (permutations > 0) {
             arg <- if (missing(strata)) NULL else strata
@@ -20,9 +38,7 @@
         permutations <- nrow(permutations)
     }
     perm <- sapply(1:permutations,
-                   function(i, ...) procrustes(X, Y[permat[i,],],
-                                               symmetric = TRUE)$ss)
-    perm <- sqrt(1 - perm)
+                   function(i, ...) procr(X, Y[permat[i,],]))
     Pval <- (sum(perm >= sol$t0) + 1)/(permutations + 1)
     if (!missing(strata)) {
         strata <- deparse(substitute(strata))

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2013-01-25 12:07:42 UTC (rev 2387)
+++ pkg/vegan/inst/ChangeLog	2013-01-25 12:51:10 UTC (rev 2388)
@@ -2,8 +2,22 @@
 
 VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
 
-Version 2.1-23 (opened January 8, 2013)
+Version 2.1-24 (opened January 25, 2013)
 
+	* protest: huge speed-up. Instead of calling procrustes() in every
+	permutation step, we only calculate the goodness of fit statistic
+	in svd(). This avoids a huge overhead of procrustes(). In a test
+	with a 160 x 12 matrix (RDA scores from Bryce Canyon data
+	'bryceveg') with 9999 permutations, the time went down from 12 sec
+	to 0.8 sec in my desktop.
+
+	* procrustes: marginally faster way of getting sum of squares of a
+	matrix. This should not influence the results, but one metaMDS()
+	trace result has a small difference in nearly-zero rmse (was
+	1.094382e-06, is 1.09439e-06) in my desktop.
+	
+Version 2.1-23 (closed January 25, 2013)
+
 	* clamtest: Richard Telford reported an extreme case 
 	{rbind(a=c(1,0,5,10,4),b=c(0,10,5,2,1))} where
 	clamtest failed due to inadequately setting up minimum