[Vegan-commits] r2887 - in pkg/vegan: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 23 09:20:05 CEST 2014
Author: jarioksa
Date: 2014-09-23 09:20:05 +0200 (Tue, 23 Sep 2014)
New Revision: 2887
Modified:
pkg/vegan/R/permutest.betadisper.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/betadisper.Rd
pkg/vegan/man/permutest.betadisper.Rd
Log:
Merge branch 'master' into r-forge-svn-local
Modified: pkg/vegan/R/permutest.betadisper.R
===================================================================
--- pkg/vegan/R/permutest.betadisper.R 2014-09-22 08:38:01 UTC (rev 2886)
+++ pkg/vegan/R/permutest.betadisper.R 2014-09-23 07:20:05 UTC (rev 2887)
@@ -1,6 +1,6 @@
`permutest.betadisper` <- function(x, pairwise = FALSE,
- permutations = how(nperm = 999),
- ...)
+ permutations = 999,
+ parallel = getOption("mc.cores"), ...)
{
t.statistic <- function(x, y) {
m <- length(x)
@@ -13,8 +13,48 @@
(xbar - ybar) / (pooled * sqrt(1/m + 1/n))
}
+ permFun <- function(idx) {
+ if (!is.matrix(idx)) {
+ dim(i) <- c(1, length(idx))
+ }
+ R <- nrow(idx)
+ Fperm <- matrix(nrow = R, ncol = 1)
+ if (pairwise) { # set up object to hold t stats
+ Tperm <- matrix(ncol = n.pairs, nrow = R)
+ Jseq <- seq_len(n.pairs)
+ }
+ rdf <- nobs - p # residual degrees of freedom
+ ## iterate
+ for (i in seq_len(R)) { # iterate
+ take <- idx[i, ] # current permutation from set
+ p.resid <- resids[take] # permute residuals
+ f <- qr.fitted(mod.Q, p.resid) # create new data
+ mss <- sum((f - mean(f))^2)
+ r <- qr.resid(mod.Q, p.resid)
+ rss <- sum(r^2)
+ resvar <- rss / rdf
+ Fperm[i, ] <- (mss / (p - 1)) / resvar
+
+ ## pairwise tests
+ if(pairwise) {
+ for(j in Jseq) {
+ grp1 <- x$distance[take][group == combin[1, j]]
+ grp2 <- x$distance[take][group == combin[2, j]]
+ Tperm[i, j] <- t.statistic(grp1, grp2)
+ }
+ }
+ }
+
+ ## bind on pairwise stats if any
+ if (pairwise) {
+ Fperm <- cbind(Fperm, Tperm)
+ }
+ Fperm
+ }
+
if(!inherits(x, "betadisper"))
stop("Only for class \"betadisper\"")
+
## will issue error if only a single group
mod.aov <- anova(x)
nobs <- length(x$distances) ## number of observations
@@ -30,52 +70,60 @@
permutations <- getPermuteMatrix(permutations, nobs)
nperm <- nrow(permutations)
- ## set-up objects to hold permuted results
- res <- numeric(length = nperm + 1)
- res[1] <- summary(mod)$fstatistic[1]
## pairwise comparisons
if(pairwise) {
- ## unique pairings
- combin <- combn(levels(x$group), 2)
+ combin <- combn(levels(x$group), 2) # unique pairings
n.pairs <- ncol(combin)
- t.stats <- matrix(0, ncol = n.pairs, nrow = nperm + 1)
- t.stats[1,] <- apply(combn(levels(group), 2), 2, function(z) {
- t.statistic(x$distances[group == z[1]],
- x$distances[group == z[2]])})
}
- ## begin loop over shuffleSet perms
- for(i in seq_len(nperm)) {
- perm <- permutations[i,] ## take current permutation from set
- perm.resid <- resids[perm] ## permute residuals
- f <- qr.fitted(mod.Q, perm.resid) ## create new data
- mss <- sum((f - mean(f))^2)
- r <- qr.resid(mod.Q, perm.resid)
- rss <- sum(r^2)
- rdf <- nobs - p
- resvar <- rss / rdf
- res[i+1] <- (mss / (p - 1)) / resvar
-
- ## pairwise comparisons
- if(pairwise) {
- for(j in seq_len(n.pairs)) {
- grp1 <- x$distance[perm][group == combin[1, j]]
- grp2 <- x$distance[perm][group == combin[2, j]]
- t.stats[i+1, j] <- t.statistic(grp1, grp2)
+ ## Parallel processing of permutations
+ if (is.null(parallel)) {
+ parallel <- 1
+ }
+ hasClus <- inherits(parallel, "cluster")
+ if ((hasClus || parallel > 1L) && requireNamespace("parallel")) {
+ if (.Platform$OS.type == "unix" && !hasClus) {
+ Pstats <- do.call("rbind",
+ mclapply(seq_len(nperm),
+ function(x) permFun(permutations[x, , drop = FALSE]),
+ mc.cores = parallel))
+ } else {
+ ## if hasClus, don't set up and top a temporary cluster
+ if (!hasClus) {
+ parallel <- makeCluster(parallel)
}
+ Pstats <- parRapply(parallel, permutations, function(x) permFun(x))
+ if (!hasClus) {
+ stopCluster(parallel)
+ }
}
+ } else {
+ Pstats <- permFun(permutations)
}
+ ## Process results
+ F0 <- summary(mod)$fstatistic[1]
+ Fstats <- round(Pstats[, 1], 12) # allow empty dim to be dropped
+ F0 <- round(F0, 12)
+ ## pairwise comparisons
+ if(pairwise) {
+ T0 <- apply(combn(levels(group), 2), 2, function(z) {
+ t.statistic(x$distances[group == z[1]],
+ x$distances[group == z[2]])})
+ Tstats <- round(Pstats[, -1, drop = FALSE], 12)
+ T0 <- round(T0, 12)
+ }
+
## compute permutation p-value
- pval <- sum(res >= res[1]) / length(res)
+ pval <- (sum(Fstats >= F0) + 1) / (length(Fstats) + 1)
if(pairwise) {
df <- apply(combin, 2, function(z) {
length(x$distances[group == z[1]]) +
length(x$distance[group == z[2]]) - 2})
- pairwise <- list(observed = 2 * pt(-abs(t.stats[1,]), df),
- permuted = apply(t.stats, 2,
- function(z) sum(abs(z) >= abs(z[1]))/length(z)))
+ pairwise <- list(observed = 2 * pt(-abs(T0), df),
+ permuted = apply(Tstats, 2,
+ function(z) (sum(abs(z) >= abs(z[1])) + 1) / (length(z) + 1)))
names(pairwise$observed) <- names(pairwise$permuted) <-
apply(combin, 2, paste, collapse = "-")
} else {
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2014-09-22 08:38:01 UTC (rev 2886)
+++ pkg/vegan/inst/ChangeLog 2014-09-23 07:20:05 UTC (rev 2887)
@@ -2,6 +2,17 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
+Version 2.1-43
+
+ * betadisper Permutation tests via the `permutest()` method
+ are now parallelised. The number of cores to parallelise over
+ is specified by argument `parallel` in common with the
+ implementation for `cca()`.
+
+ Fix a couple of bugs in the examples; number of permutations
+ was specified using `control` which is not an argument to the
+ `permutest()` method for `betadisper()`.
+
Version 2.1-42 (opened September 4, 2014)
* Opened a new version to prepare release 2.2-0.
Modified: pkg/vegan/man/betadisper.Rd
===================================================================
--- pkg/vegan/man/betadisper.Rd 2014-09-22 08:38:01 UTC (rev 2886)
+++ pkg/vegan/man/betadisper.Rd 2014-09-23 07:20:05 UTC (rev 2887)
@@ -240,7 +240,7 @@
anova(mod)
## Permutation test for F
-permutest(mod, pairwise = TRUE)
+permutest(mod, pairwise = TRUE, permutations = 99)
## Tukey's Honest Significant Differences
(mod.HSD <- TukeyHSD(mod))
@@ -279,7 +279,7 @@
dis[c(2, 20)] <- NA
mod2 <- betadisper(dis, groups) ## warnings
mod2
-permutest(mod2, control = how(nperm = 100))
+permutest(mod2, permutations = 99)
anova(mod2)
plot(mod2)
boxplot(mod2)
@@ -288,7 +288,7 @@
## Using group centroids
mod3 <- betadisper(dis, groups, type = "centroid")
mod3
-permutest(mod3, control = how(nperm = 100))
+permutest(mod3, permutations = 99)
anova(mod3)
plot(mod3)
boxplot(mod3)
Modified: pkg/vegan/man/permutest.betadisper.Rd
===================================================================
--- pkg/vegan/man/permutest.betadisper.Rd 2014-09-22 08:38:01 UTC (rev 2886)
+++ pkg/vegan/man/permutest.betadisper.Rd 2014-09-23 07:20:05 UTC (rev 2887)
@@ -10,7 +10,8 @@
}
\usage{
\method{permutest}{betadisper}(x, pairwise = FALSE,
- permutations = how(nperm = 999),
+ permutations = 999,
+ parallel = getOption("mc.cores"),
\dots)
}
%- maybe also 'usage' for other objects documented here.
@@ -22,6 +23,9 @@
as returned by the function \code{\link[permute]{how}}, or the
number of permutations required, or a permutation matrix where each
row gives the permuted indices.}
+ \item{parallel}{Number of parallel processes or a predefined socket
+ cluster. With \code{parallel = 1} uses ordinary, non-parallel
+ processing.}
\item{\dots}{Arguments passed to other methods.}
}
\details{
More information about the Vegan-commits
mailing list