[Vegan-commits] r781 - in pkg/vegan: R man

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

Author: gsimpson
Date: 2009-04-03 16:09:54 +0200 (Fri, 03 Apr 2009)
New Revision: 781

Modified:
   pkg/vegan/R/betadisper.R
   pkg/vegan/man/betadisper.Rd
Log:
Bug fix and updates to betadisper

Modified: pkg/vegan/R/betadisper.R
===================================================================
--- pkg/vegan/R/betadisper.R	2009-04-02 17:45:43 UTC (rev 780)
+++ pkg/vegan/R/betadisper.R	2009-04-03 14:09:54 UTC (rev 781)
@@ -1,9 +1,12 @@
 `betadisper` <-
-    function(d, group, type = c("centroid","median"))
+    function(d, group, type = c("centroid","median"), tol = 1e-07)
 {
     ## uses code from stats:::cmdscale by R Core Development Team
     if(!inherits(d, "dist"))
         stop("distances 'd' must be a 'dist' object")
+    type <- match.arg(type)
+    if(type == "median")
+        .NotYetUsed('type = "median"')
     ## checks for groups - need to be a factor for later
     if(!is.factor(group))
         group <- as.factor(group)
@@ -11,14 +14,23 @@
     x <- matrix(0, ncol = n, nrow = n)
     x[row(x) > col(x)] <- d^2
     x <- x + t(x)
-    .C(stats:::R_dblcen, as.double(x), as.integer(n), DUP = FALSE)
+    storage.mode(x) <- "double"
+    .C(stats:::R_dblcen, x, as.integer(n), DUP = FALSE)
     e <- eigen(-x/2, symmetric = TRUE)
     vectors <- e$vectors
     eig <- e$values
+    ## check d is Euclidean
+    w0 <- eig[n] / eig[1]
+    if(w0 > -tol)
+        r <- sum(eig > (eig[1] * tol))
+    else
+        r <- length(eig)
+    ## truncate eig if d is Euclidean
+    eig <- eig[(rs <- seq_len(r))]
+    ## scale Eigenvectors
+    vectors <- vectors[, rs, drop = FALSE] %*% diag(sqrt(abs(eig)))
     ## store which are the positive eigenvalues
     pos <- eig > 0
-    ## scale Eigenvectors
-    vectors <- vectors %*% diag(sqrt(abs(eig)))
     ## group centroids in PCoA space
     centroids <- apply(vectors, 2, function(x) tapply(x, group, mean))
     ## for each of the groups, calculate distance to centroid for
@@ -35,12 +47,12 @@
             dist.neg <- 0
         }
     } else {
-        dist.pos <- vectors[, pos, drop=FALSE] -
-            centroids[pos]
+        dist.pos <- sweep(vectors[, pos, drop=FALSE], 2,
+                          centroids[pos])
         dist.pos <- rowSums(dist.pos^2)
         if (any(!pos)) {
-            dist.neg <- vectors[, !pos, drop=FALSE] -
-                centroids[!pos]
+            dist.neg <- sweep(vectors[, !pos, drop=FALSE], 2,
+                              centroids[!pos])
             dist.neg <- rowSums(dist.neg^2)
         } else {
             dist.neg <- 0
@@ -49,7 +61,7 @@
     ## zij are the distances of each point to its group centroid
     zij <- sqrt(abs(dist.pos - dist.neg))
     ## add in correct labels
-    colnames(vectors) <- names(eig) <- paste("PCoA", 1:n, sep = "")
+    colnames(vectors) <- names(eig) <- paste("PCoA", rs, sep = "")
     if(is.matrix(centroids))
         colnames(centroids) <- names(eig)
     else

Modified: pkg/vegan/man/betadisper.Rd
===================================================================
--- pkg/vegan/man/betadisper.Rd	2009-04-02 17:45:43 UTC (rev 780)
+++ pkg/vegan/man/betadisper.Rd	2009-04-03 14:09:54 UTC (rev 781)
@@ -25,7 +25,7 @@
   Tukey's 'Honest Significant Difference' method.
 }
 \usage{
-betadisper(d, group, type = c("centroid", "median"))
+betadisper(d, group, type = c("centroid", "median"), tol = 1e-07)
 
 \method{anova}{betadisper}(object, \dots)
 
@@ -51,6 +51,7 @@
     level (i.e.~one group).}
   \item{type}{the type of analysis to perform. Only \code{type =
       "centroid"} is currently supported.}
+  \item{tol}{tolerance for checking if \code{d} is Euclidean.}
   \item{display}{character; partial match to access scores for
     \code{"sites"} or \code{"species"}.}
   \item{object, x}{an object of class \code{"betadisper"}, the result of a