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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Oct 16 19:34:50 CEST 2011


Author: jarioksa
Date: 2011-10-16 19:34:49 +0200 (Sun, 16 Oct 2011)
New Revision: 1950

Modified:
   pkg/vegan/DESCRIPTION
   pkg/vegan/R/anosim.R
   pkg/vegan/R/envfit.default.R
   pkg/vegan/R/factorfit.R
   pkg/vegan/R/mantel.R
   pkg/vegan/R/mantel.partial.R
   pkg/vegan/R/protest.R
   pkg/vegan/R/vectorfit.R
   pkg/vegan/inst/ChangeLog
   pkg/vegan/man/anosim.Rd
   pkg/vegan/man/envfit.Rd
   pkg/vegan/man/mantel.Rd
   pkg/vegan/man/procrustes.Rd
Log:
Accept matrix of 'permutations' and internally pre-adapt functions to
the 'permute' and 'parallel' packages


Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/DESCRIPTION	2011-10-16 17:34:49 UTC (rev 1950)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 2.1-3
-Date: October 4, 2011
+Version: 2.1-4
+Date: October 16, 2011
 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/anosim.R
===================================================================
--- pkg/vegan/R/anosim.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/anosim.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -30,14 +30,15 @@
     cl.vec[within] <- levels(grouping)[grouping[take]]
     cl.vec <- factor(cl.vec, levels = c("Between", levels(grouping)))
     if (permutations) {
-        perm <- rep(0, permutations)
-        for (i in 1:permutations) {
-            take <- permuted.index(N, strata)
+        ptest <- function(take, ...) {
             cl.perm <- grouping[take]
             tmp.within <- matched(irow, icol, cl.perm)
             tmp.ave <- tapply(x.rank, tmp.within, mean)
-            perm[i] <- -diff(tmp.ave)/div
+            -diff(tmp.ave)/div
         }
+        arg <- if (missing(strata)) NULL else strata
+        permat <- t(replicate(permutations, permuted.index(N, strata = arg)))
+        perm <- sapply(1:permutations, function(i) ptest(permat[i,]))
         p.val <- (1 + sum(perm >= statistic))/(1 + permutations)
         sol$signif <- p.val
         sol$perm <- perm

Modified: pkg/vegan/R/envfit.default.R
===================================================================
--- pkg/vegan/R/envfit.default.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/envfit.default.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -16,30 +16,34 @@
         env <- env[keep,, drop=FALSE]
         na.action <- structure(seq_along(keep)[!keep], class="omit")
     }
+    ## make permutation matrix for all variables handled in the next loop
+    nr <- nrow(X)
+    if (length(permutations) == 1) {
+        if (permutations > 0 ) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(nr, strata=arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != nr)
+            stop(gettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), nr))
+        permutations <- nrow(permutations)
+    }
     if (is.data.frame(env)) {
-        facts <- unlist(lapply(env, is.factor))
-        if (sum(facts)) {
+        facts <- sapply(env, is.factor)
+        if (sum(facts)) {  # have factors
             Pfac <- env[, facts, drop = FALSE]
             P <- env[, !facts, drop = FALSE]
-            if (length(P)) {
-                if (permutations) {
-                    if (!exists(".Random.seed", envir = .GlobalEnv, 
-                                inherits = FALSE)) {
-                        runif(1)
-                    }
-                    seed <- get(".Random.seed", envir = .GlobalEnv, 
-                                inherits = FALSE)
-                }
+            if (length(P)) { # also have vectors
                 vectors <- vectorfit(X, P, permutations, strata, 
                                      choices, w = w, ...)
             }
-            if (!is.null(seed)) {
-                assign(".Random.seed", seed, envir = .GlobalEnv)
-            }
             factors <- factorfit(X, Pfac, permutations, strata,
-                                 choices, w = w, ...)
+                                         choices, w = w, ...)
             sol <- list(vector = vectors, factors = factors)
-        }
+            }
         else vectors <- vectorfit(X, env, permutations, strata, 
                                   choices, w = w, ...)
     }

Modified: pkg/vegan/R/factorfit.R
===================================================================
--- pkg/vegan/R/factorfit.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/factorfit.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -1,4 +1,4 @@
-"factorfit" <-
+`factorfit` <-
     function (X, P, permutations = 0, strata, w,  ...) 
 {
     P <- as.data.frame(P)
@@ -13,16 +13,27 @@
         w <- rep(w, NR)
     r <- NULL
     pval <- NULL
-    totvar <- .C("goffactor", as.double(X), as.integer(rep(0, 
-                                                           NR)), as.double(w), as.integer(NR), as.integer(NC), as.integer(1), 
+    totvar <- .C("goffactor", as.double(X), as.integer(rep(0, NR)),
+                 as.double(w), as.integer(NR), as.integer(NC), as.integer(1), 
                  double(1), double(1), double(1), var = double(1), PACKAGE = "vegan")$var
     sol <- centroids.cca(X, P, w)
     var.id <- rep(names(P), sapply(P, nlevels))
+    ## make permutation matrix for all variables handled in the next loop
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(NR, strata=arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        permutations <- nrow(permutations)
+    }
     for (i in 1:length(P)) {
         A <- as.integer(P[[i]])
         NL <- nlevels(P[[i]])
-        invar <- .C("goffactor", as.double(X), as.integer(A - 
-                                                          1), as.double(w), as.integer(NR), as.integer(NC), 
+        invar <- .C("goffactor", as.double(X), as.integer(A - 1), as.double(w),
+                    as.integer(NR), as.integer(NC), 
                     as.integer(NL), double(NL), double(NL), double(NL), 
                     var = double(1), PACKAGE = "vegan")$var
         r.this <- 1 - invar/totvar
@@ -30,16 +41,17 @@
         if (permutations) {
             A <- as.integer(P[[i]])
             NL <- nlevels(P[[i]])
-            tmp <- rep(NA, permutations)
-            for (i in 1:permutations) {
-                indx <- permuted.index(length(A), strata)
+            ptest <- function(indx, ...) {
                 take <- A[indx]
-                invar <- .C("goffactor", as.double(X), as.integer(take - 
-                                                                  1), as.double(w), as.integer(NR), as.integer(NC), 
+                invar <- .C("goffactor", as.double(X),
+                            as.integer(take -  1), as.double(w),
+                            as.integer(NR), as.integer(NC), 
                             as.integer(NL), double(NL), double(NL), double(NL), 
                             var = double(1), PACKAGE = "vegan")$var
-                tmp[i] <- 1 - invar/totvar
+                1 - invar/totvar
             }
+            tmp <- sapply(1:permutations,
+                          function(indx,...) ptest(permat[indx,], ...))
             pval.this <- (sum(tmp > r.this) + 1)/(permutations + 1)
             pval <- c(pval, pval.this)
         }

Modified: pkg/vegan/R/mantel.R
===================================================================
--- pkg/vegan/R/mantel.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/mantel.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -7,18 +7,31 @@
     tmp <- cor.test(as.vector(xdis), ydis, method = method)
     statistic <- as.numeric(tmp$estimate)
     variant <- tmp$method
+    N <- attr(xdis, "Size")
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(N, strata = arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != N)
+            stop(gettextf("'permutations' have %d columns, but data have %d observations",
+                          ncol(permat), N))
+        permutations <- nrow(permutations)
+    }
     if (permutations) {
-        N <- attr(xdis, "Size")
-        perm <- rep(0, permutations)
-        ## asdist asn an index selects lower diagonal like as.dist,
-        ## but is faster since it does not seet 'dist' attributes
+        perm <- numeric(permutations)
+        ## asdist as an index selects lower diagonal like as.dist,
+        ## but is faster since it does not set 'dist' attributes
         xmat <- as.matrix(xdis)
         asdist <- row(xmat) > col(xmat)
-        for (i in 1:permutations) {
-            take <- permuted.index(N, strata)
+        ptest <- function(take, ...) {
             permvec <- (xmat[take, take])[asdist]
-            perm[i] <- cor(permvec, ydis, method = method)
+            drop(cor(permvec, ydis, method = method))
         }
+        perm <- sapply(1:permutations, function(i, ...) ptest(permat[i,], ...) )
         signif <- (sum(perm >= statistic) + 1)/(permutations + 1)
      }
     else {

Modified: pkg/vegan/R/mantel.partial.R
===================================================================
--- pkg/vegan/R/mantel.partial.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/mantel.partial.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -14,18 +14,32 @@
     variant <- rxy$method
     rxy <- rxy$estimate
     statistic <- part.cor(rxy, rxz, ryz)
+    N <- attr(xdis, "Size")
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if(missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(N, strata = arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != N)
+            stop(gettextf("'permutations' have %d columns, but data have %d observations",
+                          ncol(permat), N))
+        permutations <- nrow(permutations)
+    }
     if (permutations) {
         N <- attr(xdis, "Size")
         perm <- rep(0, permutations)
         xmat <- as.matrix(xdis)
         asdist <- row(xmat) > col(xmat)
-        for (i in 1:permutations) {
-            take <- permuted.index(N, strata)
+        ptest <- function(take, ...) {
             permvec <- (xmat[take, take])[asdist]
             rxy <- cor(permvec, ydis, method = method)
             rxz <- cor(permvec, zdis, method = method)
-            perm[i] <- part.cor(rxy, rxz, ryz)
+            part.cor(rxy, rxz, ryz)
         }
+        perm <- sapply(1:permutations, function(i, ...) ptest(permat[i,], ...))
         signif <- (sum(perm >= statistic)+1)/(permutations + 1)
     }
     else {

Modified: pkg/vegan/R/protest.R
===================================================================
--- pkg/vegan/R/protest.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/protest.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -1,4 +1,4 @@
-"protest" <-
+`protest` <-
     function (X, Y, scores = "sites", permutations = 999, strata, ...)
 {
     X <- scores(X, display = scores, ...)
@@ -6,15 +6,25 @@
     sol <- procrustes(X, Y, symmetric = TRUE)
     sol$t0 <- sqrt(1 - sol$ss)
     N <- nrow(X)
-    perm <- rep(0, permutations)
-    for (i in 1:permutations) {
-        take <- permuted.index(N, strata)
-        tmp <- procrustes(X, Y[take, ], symmetric = TRUE)$ss
-        perm[i] <- sqrt(1 - tmp)
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(N, strata = arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != N)
+            stop(gettextf("'permutations' have %d columns, but data have %d observations",
+                          ncol(permat), N))
+        permutations <- nrow(permutations)
     }
+    perm <- sapply(1:permutations,
+                   function(i, ...) procrustes(X, Y[permat[i,],],
+                                               symmetric = TRUE)$ss)
+    perm <- sqrt(1 - perm)
     perm <- c(sol$t0, perm)
-    permutations <- permutations + 1
-    Pval <- sum(perm >= sol$t0)/permutations
+    Pval <- sum(perm >= sol$t0)/(permutations + 1)
     if (!missing(strata)) {
         strata <- deparse(substitute(strata))
         s.val <- strata

Modified: pkg/vegan/R/vectorfit.R
===================================================================
--- pkg/vegan/R/vectorfit.R	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/R/vectorfit.R	2011-10-16 17:34:49 UTC (rev 1950)
@@ -23,21 +23,38 @@
     if (is.null(colnames(X))) 
         colnames(heads) <- paste("Dim", 1:nc, sep = "")
     else colnames(heads) <- colnames(X)
+    ## make permutation matrix for all variables handled in the next loop
+    nr <- nrow(X)
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if(missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(nr, strata = arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != nr)
+            stop(gettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), nr))
+        permutations <- nrow(permutations)
+    }
     if (permutations) {
-        nr <- nrow(X)
-        permstore <- matrix(nrow = permutations, ncol = ncol(P))
-        for (i in 1:permutations) {
-            indx <- permuted.index(nrow(P), strata)
+        ptest <- function(indx, ...) {
             take <- P[indx, , drop = FALSE]
             take <- .C("wcentre", x = as.double(take), as.double(w),
                        as.integer(nrow(take)), as.integer(ncol(take)),
                        PACKAGE = "vegan")$x
             dim(take) <- dim(P)
             Hperm <- qr.fitted(Q, take)
-            permstore[i, ] <- diag(cor(Hperm, take))^2
+            diag(cor(Hperm, take))^2
         }
-        permstore <- sweep(permstore, 2, r, ">")
-        pvals <- (apply(permstore, 2, sum) + 1)/(permutations + 1)
+        permstore <- sapply(1:permutations, function(indx, ...) ptest(permat[indx,], ...))
+        ## Single variable is dropped to a vector, and otherwise
+        ## permutations are the matrix columns and variables are rows
+        if (!is.matrix(permstore))
+            permstore <- matrix(permstore, ncol=permutations)
+        permstore <- sweep(permstore, 1, r, ">")
+        pvals <- (rowSums(permstore) + 1)/(permutations + 1)
     }
     else pvals <- NULL
     sol <- list(arrows = heads, r = r, permutations = permutations, 

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/inst/ChangeLog	2011-10-16 17:34:49 UTC (rev 1950)
@@ -2,8 +2,23 @@
 
 VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
 
-Version 2.1-3 (opened October 4, 2011)
+Version 2.1-4 (opened October 16, 2011)
 
+	* anosim, envfit (factorfit, vectorfit), mantel, mantel.partial,
+	protest: user interface changed and 'permutations' can now be a
+	matrix where each row gives permuted indices. Internally first
+	find a permutation matrix or use the given permutation matrix, and
+	then find the statistics with single {ls}apply. This makes the
+	functions ready both for the 'permute' package and for
+	parallelization (replace {ls}apply with mclapply, par{SL}apply).
+	Function envfit() was much simplified by generating a common
+	permutation matrix in envfit.default() and using that as the input
+	to vectorfit() and factorfit(). Functions CCorA not yet
+	changed. The anova.cca* cases should also be made to use a single
+	generated permutation matrix, as permutest.cca() allows this.
+
+Version 2.1-3 (closed October 16, 2011)
+
 	* added plot.preston, lines.preston and plot.fisher (that also can
 	add points and lines). These are similar as plot.prestonfit and
 	plot.fisherfit, but without the fitted model. Among other things,

Modified: pkg/vegan/man/anosim.Rd
===================================================================
--- pkg/vegan/man/anosim.Rd	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/man/anosim.Rd	2011-10-16 17:34:49 UTC (rev 1950)
@@ -18,8 +18,8 @@
     columns are response variable(s), or a dissimilarity object or a
     symmetric square matrix of dissimilarities.}
   \item{grouping}{Factor for grouping observations.}
-  \item{permutations}{Number of permutation to assess the significance
-    of the ANOSIM statistic. }
+  \item{permutations}{Number of permutations or a permutation matrix
+    where each row gives the permuted indices.}
   \item{distance}{Choice of distance metric that measures the
     dissimilarity between two observations . See \code{\link{vegdist}} for
     options.  This will be used if \code{dat} was not a dissimilarity

Modified: pkg/vegan/man/envfit.Rd
===================================================================
--- pkg/vegan/man/envfit.Rd	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/man/envfit.Rd	2011-10-16 17:34:49 UTC (rev 1950)
@@ -36,8 +36,9 @@
   \item{P}{Data frame, matrix or vector of environmental
     variable(s). These must be continuous for \code{vectorfit} and
     factors or characters for \code{factorfit}. }
-  \item{permutations}{ Number of permutations for assessing significance
-    of vectors or factors. Set to \code{0} to skip permutations.}
+  \item{permutations}{Number of permutations or a permutation matrix
+    where each row gives the permuted indices. Set \code{permutations = 0}
+    to skip the permutations.}
   \item{formula, data}{Model  \code{\link{formula}} and data.  }
   \item{na.rm}{Remove points with missing values in ordination scores
     or environmental variables. The operation is casewise: the whole

Modified: pkg/vegan/man/mantel.Rd
===================================================================
--- pkg/vegan/man/mantel.Rd	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/man/mantel.Rd	2011-10-16 17:34:49 UTC (rev 1950)
@@ -22,7 +22,8 @@
   \item{xdis, ydis, zdis}{ Dissimilarity matrices or a \code{dist} objects. }
   \item{method}{ Correlation method, as accepted by \code{\link{cor}}:
     \code{"pearson"}, \code{"spearman"} or \code{"kendall"}. }
-  \item{permutations}{Number of permutations in assessing significance. }
+  \item{permutations}{Number of permutations or a permutation matrix
+    where each row gives the permuted indices.}
   \item{strata}{An integer vector or factor specifying the strata for
     permutation. If supplied, observations are permuted only within the
     specified strata.}

Modified: pkg/vegan/man/procrustes.Rd
===================================================================
--- pkg/vegan/man/procrustes.Rd	2011-10-11 19:00:46 UTC (rev 1949)
+++ pkg/vegan/man/procrustes.Rd	2011-10-16 17:34:49 UTC (rev 1950)
@@ -65,8 +65,9 @@
     \code{truemean = FALSE}.}
   \item{newdata}{Matrix of coordinates to be rotated and translated to
      the target.}
-  \item{permutations}{Number of permutation to assess the significance
-    of the symmetric Procrustes statistic. }
+  \item{permutations}{Number of permutations or a permutation matrix
+    where each row gives the permuted indices. These are used to asses the
+    signficance of the symmetric Procrustes statistic.}
   \item{strata}{An integer vector or factor specifying the strata for
     permutation. If supplied, observations are permuted only within the
     specified strata.}



More information about the Vegan-commits mailing list