[Vegan-commits] r2602 - pkg/vegan/R

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

Author: jarioksa
Date: 2013-09-05 12:00:25 +0200 (Thu, 05 Sep 2013)
New Revision: 2602

Modified:
   pkg/vegan/R/MDSrotate.R
Log:
MDSrotate handles more than one external vector

Modified: pkg/vegan/R/MDSrotate.R
===================================================================
--- pkg/vegan/R/MDSrotate.R	2013-09-05 09:57:29 UTC (rev 2601)
+++ pkg/vegan/R/MDSrotate.R	2013-09-05 10:00:25 UTC (rev 2602)
@@ -15,46 +15,43 @@
     N <- NCOL(x)
     if (N < 2)
         stop(gettextf("needs at least 2 dimensions"))
-    vec <- drop(vec)
-    if (length(dim(vec)) > 1)
-        stop(gettextf("function works only with univariate 'vec'"))
+    vec <- as.matrix(vec)
+    NV <- NCOL(vec)
+    if (NV >= N)
+        stop(gettextf("You can have max %d vectors, but you had %d",
+             N-1, NV))
     if (!is.numeric(vec))
         stop(gettextf("'vec' must be numeric"))
     ## vectorfit finds the direction cosine. We rotate first axis to
     ## 'vec' which means that we make other axes orthogonal to 'vec'
     ## one by one
     if (na.rm)
-        keep <- !is.na(vec)
+        keep <- complete.cases(vec)
     else
-        keep <- !logical(length(vec))
-    ## scores must be orthogonal for the next loop to work
-    if (N > 2) {
-        pc <- prcomp(x[keep,])
-        x <- x %*% pc$rotation
-        if (!all(is.na(sp)))
-            sp <- sp %*% pc$rotation
-    }
+        keep <- !logical(NROW(vec))
     ## Rotation loop
-    for (k in 2:N) {
-        arrs <- vectorfit(x[keep,], vec[keep], permutations = 0)$arrows
-        rot <- arrs[c(1,k)]/sqrt(sum(arrs[c(1,k)]^2))
-        rot <- drop(rot)
-        ## counterclockwise rotation matrix:
-        ## [cos theta   -sin theta]
-        ## [sin theta    cos theta]
-        rot <- rbind(rot, rev(rot))
-        rot[1,2] <- -rot[1,2]
-        ## Rotation of points and species scores
-        x[, c(1,k)] <- x[, c(1,k)] %*% rot
-        if (!all(is.na(sp)))
-            sp[, c(1,k)] <- sp[, c(1,k)] %*% rot
+    for(v in seq_len(NV)) {
+        for (k in (v+1):N) {
+            arrs <- vectorfit(x[keep,], vec[keep,v], permutations = 0)$arrows
+            rot <- arrs[c(v,k)]/sqrt(sum(arrs[c(v,k)]^2))
+            rot <- drop(rot)
+            ## counterclockwise rotation matrix:
+            ## [cos theta   -sin theta]
+            ## [sin theta    cos theta]
+            rot <- rbind(rot, rev(rot))
+            rot[1,2] <- -rot[1,2]
+            ## Rotation of points and species scores
+            x[, c(v,k)] <- x[, c(v,k)] %*% rot
+            if (!all(is.na(sp)))
+                sp[, c(v,k)] <- sp[, c(v,k)] %*% rot
+        }
     }
-    ## Rotate 2..N axes to PC
-    if (N > 2 && attr(object$points, "pc")) {
-        pc <- prcomp(x[,-1])
-        x[,-1] <- pc$x
+    ## Two or more free axes are (optionally) rotated to PCs
+    if (N - NV > 1 && attr(object$points, "pc")) {
+        pc <- prcomp(x[,-seq_len(NV)])
+        x[,-seq_len(NV)] <- pc$x
         if (!all(is.na(sp)))
-            sp[,-1] <- sp[,-1] %*% pc$rotation
+            sp[,-seq_len(NV)] <- sp[,-seq_len(NV)] %*% pc$rotation
     }
     ## '[] <-' retains attributes
     object$points[] <- x