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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Oct 18 15:37:37 CEST 2011


Author: jarioksa
Date: 2011-10-18 15:37:36 +0200 (Tue, 18 Oct 2011)
New Revision: 1953

Modified:
   pkg/vegan/R/CCorA.R
   pkg/vegan/R/adonis.R
   pkg/vegan/R/anosim.R
   pkg/vegan/R/mrpp.R
   pkg/vegan/inst/ChangeLog
   pkg/vegan/man/CCorA.Rd
   pkg/vegan/man/adonis.Rd
   pkg/vegan/man/mrpp.Rd
Log:
CCorA, adonis and mrpp can take a permutation matrix as input,
and CCorA was made 'permute'/'parallel' ready (mrpp, adonis already were)


Modified: pkg/vegan/R/CCorA.R
===================================================================
--- pkg/vegan/R/CCorA.R	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/R/CCorA.R	2011-10-18 13:37:36 UTC (rev 1953)
@@ -40,20 +40,15 @@
         }
         invisible(0)
     }
-    probPillai <- function(Y, X, n, S11.inv, S22.inv, s, df1, df2, epsilon,
-                           Fref, nperm, ...) {
+    probPillai <- function(Y.per, X, n, S11.inv, S22.inv, s, df1, df2, epsilon,
+                           Fref, permat, ...) {
         ## Permutation test for Pillai's trace in CCorA.
         ## Reference: Brian McArdle's unpublished graduate course notes.
-        nGE <- 1
-        for(i in 1:nperm) {
-            Y.per <- Y[permuted.index(n, ...),, drop=FALSE]
-            S12.per <- cov(Y.per,X)
-            gross.mat <- S12.per %*% S22.inv %*% t(S12.per) %*% S11.inv
-            Pillai.per <- sum(diag(gross.mat))
-            Fper  <- (Pillai.per*df2)/((s-Pillai.per)*df1)
-            if(Fper >= (Fref-epsilon)) nGE <- nGE+1
-        }
-        P <- nGE/(nperm+1)
+        S12.per <- cov(Y.per,X)
+        gross.mat <- S12.per %*% S22.inv %*% t(S12.per) %*% S11.inv
+        Pillai.per <- sum(diag(gross.mat))
+        Fper  <- (Pillai.per*df2)/((s-Pillai.per)*df1)
+        Fper >= (Fref-epsilon)
     }
     ## END: internal functions
     ##
@@ -161,14 +156,25 @@
     df2 <- (n - max(pp,qq) - 1)
     Fval  <- (PillaiTrace*df2)/((s-PillaiTrace)*df1)
     p.Pillai <- pf(Fval, s*df1, s*df2, lower.tail=FALSE)
-
-    if(nperm > 0) {
-        p.perm <- probPillai(Y, X, n, S11.inv, S22.inv, s, df1, df2,
-                             epsilon, Fval, nperm, ...)
+    if (length(nperm) == 1) {
+        if (nperm > 0)
+            permat <- t(replicate(nperm, permuted.index(n, ...)))
+    } else  {
+        permat <- as.matrix(nperm)
+        if (ncol(permat) != n)
+            stop(gettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), n))
+        nperm <- nrow(permat)
+    }
+    if (nperm > 0) {
+        p.perm <- sapply(1:nperm, function(indx, ...) 
+                         probPillai(Y[permat[indx,],] , X, n, S11.inv, S22.inv, s,
+                                    df1, df2, epsilon, Fval, nperm, ...))
+        p.perm <- (sum(p.perm) +1)/(nperm + 1)
     } else {
         p.perm <- NA
     }
-
+    
     out <- list(Pillai=PillaiTrace, Eigenvalues=Eigenvalues, CanCorr=K.svd$d,
                 Mat.ranks=c(RsquareX.Y$m, RsquareY.X$m), 
                 RDA.Rsquares=c(RsquareY.X$Rsquare, RsquareX.Y$Rsquare),

Modified: pkg/vegan/R/adonis.R
===================================================================
--- pkg/vegan/R/adonis.R	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/R/adonis.R	2011-10-18 13:37:36 UTC (rev 1953)
@@ -91,11 +91,18 @@
           ) }
     
     ## Permutations
-    if (missing(strata)) 
-        strata <- NULL
-    p <- sapply(1:permutations,
-                function(x) permuted.index(n, strata=strata))
-
+    if (length(permutations) == 1) {
+        if (missing(strata)) 
+            strata <- NULL
+        p <- replicate(permutations,
+                       permuted.index(n, strata=strata))
+    } else {
+        p <- t(as.matrix(permutations))
+        if (nrow(p) != n)
+            stop(gettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), n))
+        permutations <- ncol(p)
+    }
     
     tH.s <- sapply(H.s, t)
     tIH.snterm <- t(I-H.snterm)

Modified: pkg/vegan/R/anosim.R
===================================================================
--- pkg/vegan/R/anosim.R	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/R/anosim.R	2011-10-18 13:37:36 UTC (rev 1953)
@@ -29,20 +29,28 @@
     take <- as.numeric(irow[within])
     cl.vec[within] <- levels(grouping)[grouping[take]]
     cl.vec <- factor(cl.vec, levels = c("Between", levels(grouping)))
-    if (permutations) {
-        ptest <- function(take, ...) {
-            cl.perm <- grouping[take]
-            tmp.within <- matched(irow, icol, cl.perm)
-            tmp.ave <- tapply(x.rank, tmp.within, mean)
-            -diff(tmp.ave)/div
+    ptest <- function(take, ...) {
+        cl.perm <- grouping[take]
+        tmp.within <- matched(irow, icol, cl.perm)
+        tmp.ave <- tapply(x.rank, tmp.within, mean)
+        -diff(tmp.ave)/div
+    }
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations, permuted.index(N, strata = arg)))
         }
-        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
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != N)
+            stop(gettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), N))
+        permutations <- nrow(permat)
     }
+    perm <- sapply(1:permutations, function(i) ptest(permat[i,]))
+    p.val <- (1 + sum(perm >= statistic))/(1 + permutations)
+    sol$signif <- p.val
+    sol$perm <- perm
     sol$permutations <- permutations
     sol$statistic <- as.numeric(statistic)
     sol$class.vec <- cl.vec

Modified: pkg/vegan/R/mrpp.R
===================================================================
--- pkg/vegan/R/mrpp.R	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/R/mrpp.R	2011-10-18 13:37:36 UTC (rev 1953)
@@ -37,11 +37,18 @@
     ## significance test for it. Keep the item in reserve for
     ## possible later re-inclusion.
     CS <- NA
-    if (missing(strata)) 
-        strata <- NULL
-    perms <- sapply(1:permutations, function(x) grouping[permuted.index(N, 
-        strata = strata)])
-    m.ds <- numeric(permutations)
+    if (length(permutations) == 1) {
+        if (missing(strata)) 
+            strata <- NULL
+        perms <- sapply(1:permutations,
+                        function(x) grouping[permuted.index(N, strata = strata)])
+    } else {
+        perms <- apply(permutations, 1, function(indx) grouping[indx])
+        permutations <- ncol(perms)
+        if (nrow(perms) != N)
+            stop(fettextf("'permutations' have %d columns, but data have %d rows",
+                          ncol(permat), n))
+    }
     m.ds <- apply(perms, 2, function(x) mrpp.perms(x, dmat, indls, 
         w))
     p <- (1 + sum(del >= m.ds))/(permutations + 1)
@@ -50,7 +57,7 @@
         n = ncl, classdelta = classdel,
         Pvalue = p, A = r2, distance = distance, weight.type = weight.type, 
         boot.deltas = m.ds, permutations = permutations)
-    if (!is.null(strata)) {
+    if (!missing(strata) && !is.null(strata)) {
         out$strata <- deparse(substitute(strata))
         out$stratum.values <- strata
     }

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/inst/ChangeLog	2011-10-18 13:37:36 UTC (rev 1953)
@@ -4,18 +4,20 @@
 
 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.
+	* adonis, anosim, CCorA, envfit (factorfit, vectorfit), mantel,
+	mantel.partial, mrpp, 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. Functions adonis and mrpp already worked like this, but
+	they gained the option of matrix input. 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(). 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)
 

Modified: pkg/vegan/man/CCorA.Rd
===================================================================
--- pkg/vegan/man/CCorA.Rd	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/man/CCorA.Rd	2011-10-18 13:37:36 UTC (rev 1953)
@@ -23,7 +23,8 @@
   \item{stand.Y}{ Logical; should \code{Y} be standardized? }
   \item{stand.X}{ Logical; should \code{X} be standardized? }
   \item{nperm}{ Numeric; number of permutations to evaluate the
-    significance of Pillai's trace, e.g. \code{nperm=99} or \code{nperm=999}.}
+    significance of Pillai's trace, or a permutation matrix
+    where each row gives the permuted indices.} 
   \item{x}{\code{CCoaR} result object.}
   \item{plot.type}{ A character string indicating which of the following 
     plots should be produced: \code{"objects"}, \code{"variables"}, \code{"ov"} 

Modified: pkg/vegan/man/adonis.Rd
===================================================================
--- pkg/vegan/man/adonis.Rd	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/man/adonis.Rd	2011-10-18 13:37:36 UTC (rev 1953)
@@ -24,8 +24,9 @@
   Value below).} 
   \item{data}{ the data frame from which \code{A}, \code{B}, and
     \code{C} would be drawn.} 
-  \item{permutations}{ number of replicate permutations used for the
-    hypothesis tests (\eqn{F} tests).} 
+  \item{permutations}{ number of replicate permutations or a
+    permutation matrix where each row gives the permuted indices used
+    for the hypothesis tests (\eqn{F} tests).}
   \item{method}{ the name of any method used in \code{\link{vegdist}} to
     calculate pairwise distances if the left hand side of the
     \code{formula} was a data frame or a matrix. } 

Modified: pkg/vegan/man/mrpp.Rd
===================================================================
--- pkg/vegan/man/mrpp.Rd	2011-10-18 07:15:22 UTC (rev 1952)
+++ pkg/vegan/man/mrpp.Rd	2011-10-18 13:37:36 UTC (rev 1953)
@@ -25,8 +25,9 @@
     columns are response variable(s), or a dissimilarity object or a
     symmetric square matrix of dissimilarities.} 
   \item{grouping}{ Factor or numeric index for grouping observations.}
-  \item{permutations}{Number of permutations to assess the significance
-    of the MRPP statistic, \eqn{delta}.} 
+  \item{permutations}{Number of permutations or a permutation matrix
+    where each row gives the permuted indices. These are used to assess 
+    the significance of the MRPP statistic, \eqn{delta}.} 
   \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



More information about the Vegan-commits mailing list