[Vegan-commits] r2957 - in pkg/vegan: R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 28 10:22:19 CEST 2015
Author: jarioksa
Date: 2015-08-28 10:22:19 +0200 (Fri, 28 Aug 2015)
New Revision: 2957
Modified:
pkg/vegan/R/adonis.R
pkg/vegan/R/anosim.R
pkg/vegan/R/anova.cca.R
pkg/vegan/R/anova.ccabyterm.R
pkg/vegan/R/anova.ccalist.R
pkg/vegan/R/factorfit.R
pkg/vegan/R/mantel.R
pkg/vegan/R/mantel.partial.R
pkg/vegan/R/mrpp.R
pkg/vegan/R/mso.R
pkg/vegan/R/oecosimu.R
pkg/vegan/R/ordiareatest.R
pkg/vegan/R/permutest.betadisper.R
pkg/vegan/R/permutest.cca.R
pkg/vegan/R/print.permutest.cca.R
pkg/vegan/R/protest.R
pkg/vegan/R/simper.R
pkg/vegan/R/vectorfit.R
pkg/vegan/inst/NEWS.Rd
Log:
Merge branch 'cran-2.3' into r-forge-svn-local
Modified: pkg/vegan/R/adonis.R
===================================================================
--- pkg/vegan/R/adonis.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/adonis.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -3,6 +3,7 @@
contr.unordered="contr.sum", contr.ordered="contr.poly",
parallel = getOption("mc.cores"), ...)
{
+ EPS <- sqrt(.Machine$double.eps) ## use with >= in permutation P-values
## formula is model formula such as Y ~ A + B*C where Y is a data
## frame or a matrix, and A, B, and C may be factors or continuous
## variables. data is the data frame from which A, B, and C would
@@ -132,10 +133,7 @@
## Close socket cluster if created here
if (isParal && !isMulticore && !hasClus)
stopCluster(parallel)
- ## Round to avoid arbitrary P-values with tied data
- f.perms <- round(f.perms, 12)
- F.Mod <- round(F.Mod, 12)
- P <- (rowSums(t(f.perms) >= F.Mod)+1)/(permutations+1)
+ P <- (rowSums(t(f.perms) >= F.Mod - EPS)+1)/(permutations+1)
} else { # no permutations
f.perms <- P <- rep(NA, nterms)
}
Modified: pkg/vegan/R/anosim.R
===================================================================
--- pkg/vegan/R/anosim.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/anosim.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function (dat, grouping, permutations = 999,
distance = "bray", strata = NULL, parallel = getOption("mc.cores"))
{
+ EPS <- sqrt(.Machine$double.eps)
if (inherits(dat, "dist"))
x <- dat
else if (is.matrix(dat) && nrow(dat) == ncol(dat) && all(dat[lower.tri(dat)] ==
@@ -65,7 +66,7 @@
} else {
perm <- sapply(1:permutations, function(i) ptest(permat[i,]))
}
- p.val <- (1 + sum(perm >= statistic))/(1 + permutations)
+ p.val <- (1 + sum(perm >= statistic - EPS))/(1 + permutations)
} else { # no permutations
p.val <- perm <- NA
}
Modified: pkg/vegan/R/anova.cca.R
===================================================================
--- pkg/vegan/R/anova.cca.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/anova.cca.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -4,6 +4,7 @@
parallel = getOption("mc.cores"), strata = NULL,
cutoff = 1, scope = NULL)
{
+ EPS <- sqrt(.Machine$double.eps) # for permutation P-values
model <- match.arg(model)
## permutation matrix
N <- nrow(object$CA$u)
@@ -58,7 +59,7 @@
tst <- permutest.cca(object, permutations = permutations,
model = model, parallel = parallel, ...)
Fval <- c(tst$F.0, NA)
- Pval <- (sum(tst$F.perm >= tst$F.0) + 1)/(tst$nperm + 1)
+ Pval <- (sum(tst$F.perm >= tst$F.0 - EPS) + 1)/(tst$nperm + 1)
Pval <- c(Pval, NA)
table <- data.frame(tst$df, tst$chi, Fval, Pval)
if (inherits(object, "capscale") &&
Modified: pkg/vegan/R/anova.ccabyterm.R
===================================================================
--- pkg/vegan/R/anova.ccabyterm.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/anova.ccabyterm.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -56,6 +56,7 @@
`anova.ccabymargin` <-
function(object, permutations, scope, ...)
{
+ EPS <- sqrt(.Machine$double.eps)
nperm <- nrow(permutations)
## Refuse to handle models with missing data
if (!is.null(object$na.action))
@@ -95,7 +96,7 @@
Fval <- sweep(Fval, 2, Df, "/")
Fval <- sweep(Fval, 1, scale, "/")
## Simulated P-values
- Pval <- (colSums(sweep(Fval, 2, Fstat, ">=")) + 1)/(nperm + 1)
+ Pval <- (colSums(sweep(Fval, 2, Fstat - EPS, ">=")) + 1)/(nperm + 1)
## Collect results to anova data.frame
out <- data.frame(c(Df, dfbig), c(Chisq, chibig),
c(Fstat, NA), c(Pval, NA))
@@ -123,6 +124,7 @@
`anova.ccabyaxis` <-
function(object, permutations, model, parallel, cutoff = 1)
{
+ EPS <- sqrt(.Machine$double.eps)
nperm <- nrow(permutations)
## Observed F-values and Df
eig <- object$CCA$eig
Modified: pkg/vegan/R/anova.ccalist.R
===================================================================
--- pkg/vegan/R/anova.ccalist.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/anova.ccalist.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -1,6 +1,7 @@
`anova.ccalist` <-
function(object, permutations, model, parallel)
{
+ EPS <- sqrt(.Machine$double.eps)
## 'object' *must* be a list of cca objects, and 'permutations'
## *must* be a permutation matrix -- we assume that calling
## function takes care of this, and this function is not directly
@@ -67,8 +68,8 @@
pfvals <- matrix(pfvals, nrow=1, ncol=nperm)
pfvals <- sweep(pfvals, 1, df, "/")
pfvals <- sweep(pfvals, 2, pscale, "/")
- pval <- rowSums(sweep(pfvals, 1, fval, ">="))
- pval <- (pval + 1)/(nperm+1)
+ pval <- rowSums(sweep(pfvals, 1, fval - EPS, ">="))
+ pval <- (pval + 1)/(nperm + 1)
## collect table
table <- data.frame(resdf, resdev, c(NA, df),
c(NA,changedev), c(NA,fval), c(NA,pval))
Modified: pkg/vegan/R/factorfit.R
===================================================================
--- pkg/vegan/R/factorfit.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/factorfit.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -1,6 +1,7 @@
`factorfit` <-
function (X, P, permutations = 0, strata = NULL, w, ...)
{
+ EPS <- sqrt(.Machine$double.eps)
P <- as.data.frame(P)
## Check that all variables are factors, and coerce if necessary
if(any(!sapply(P, is.factor)))
@@ -50,7 +51,7 @@
}
tmp <- sapply(seq_len(permutations),
function(indx,...) ptest(permat[indx,], ...))
- pval.this <- (sum(tmp >= r.this) + 1)/(permutations + 1)
+ pval.this <- (sum(tmp >= r.this - EPS) + 1)/(permutations + 1)
pval <- c(pval, pval.this)
}
}
Modified: pkg/vegan/R/mantel.R
===================================================================
--- pkg/vegan/R/mantel.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/mantel.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function (xdis, ydis, method = "pearson", permutations = 999,
strata = NULL, na.rm = FALSE, parallel = getOption("mc.cores"))
{
+ EPS <- sqrt(.Machine$double.eps)
xdis <- as.dist(xdis)
ydis <- as.vector(as.dist(ydis))
## Handle missing values
@@ -54,7 +55,7 @@
} else {
perm <- sapply(1:permutations, function(i, ...) ptest(permat[i,], ...))
}
- signif <- (sum(perm >= statistic) + 1)/(permutations + 1)
+ signif <- (sum(perm >= statistic - EPS) + 1)/(permutations + 1)
}
else {
signif <- NA
Modified: pkg/vegan/R/mantel.partial.R
===================================================================
--- pkg/vegan/R/mantel.partial.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/mantel.partial.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function (xdis, ydis, zdis, method = "pearson", permutations = 999,
strata = NULL, na.rm = FALSE, parallel = getOption("mc.cores"))
{
+ EPS <- sqrt(.Machine$double.eps)
part.cor <- function(rxy, rxz, ryz) {
(rxy - rxz * ryz)/sqrt(1-rxz*rxz)/sqrt(1-ryz*ryz)
}
@@ -62,7 +63,7 @@
} else {
perm <- sapply(1:permutations, function(i, ...) ptest(permat[i,], ...))
}
- signif <- (sum(perm >= statistic)+1)/(permutations + 1)
+ signif <- (sum(perm >= statistic - EPS)+1)/(permutations + 1)
}
else {
signif <- NA
Modified: pkg/vegan/R/mrpp.R
===================================================================
--- pkg/vegan/R/mrpp.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/mrpp.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -3,6 +3,7 @@
weight.type = 1, strata = NULL,
parallel = getOption("mc.cores"))
{
+ EPS <- sqrt(.Machine$double.eps)
classmean <- function(ind, dmat, indls) {
sapply(indls, function(x)
mean(c(dmat[ind == x, ind == x]),
@@ -69,7 +70,7 @@
} else {
m.ds <- apply(perms, 2, function(x) mrpp.perms(x, dmat, indls, w))
}
- p <- (1 + sum(del >= m.ds))/(permutations + 1)
+ p <- (1 + sum(del + EPS >= m.ds))/(permutations + 1)
r2 <- 1 - del/E.del
} else { # no permutations
m.ds <- p <- r2 <- NA
Modified: pkg/vegan/R/mso.R
===================================================================
--- pkg/vegan/R/mso.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/mso.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function (object.cca, object.xy, grain = 1, round.up = FALSE,
permutations = 0)
{
+ EPS <- sqrt(.Machine$double.eps)
if (inherits(object.cca, "mso")) {
rm <- which(class(object.cca) == "mso")
class(object.cca) <- class(object.cca)[-rm]
@@ -76,7 +77,8 @@
}
perm <- sapply(1:nperm, function(take) permfunc(permat[take,]))
object$vario$CA.signif <-
- (rowSums(sweep(perm, 1, statistic, ">=")) + 1)/(nperm + 1)
+ (rowSums(sweep(perm, 1, statistic - EPS, ">=")) + 1)/
+ (nperm + 1)
attr(object$vario, "control") <- attr(permat, "control")
}
object$call <- match.call()
Modified: pkg/vegan/R/oecosimu.R
===================================================================
--- pkg/vegan/R/oecosimu.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/oecosimu.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -143,8 +143,15 @@
z <- (indstat - means)/sd
if (any(sd < sqrt(.Machine$double.eps)))
z[sd < sqrt(.Machine$double.eps)] <- 0
- pless <- rowSums(indstat >= simind, na.rm = TRUE)
- pmore <- rowSums(indstat <= simind, na.rm = TRUE)
+ ## results can be integers or real: comparisons differ
+ if (is.integer(indstat) && is.integer(simind)) {
+ pless <- rowSums(indstat >= simind, na.rm = TRUE)
+ pmore <- rowSums(indstat <= simind, na.rm = TRUE)
+ } else {
+ EPS <- sqrt(.Machine$double.eps)
+ pless <- rowSums(indstat + EPS >= simind, na.rm = TRUE)
+ pmore <- rowSums(indstat - EPS <= simind, na.rm = TRUE)
+ }
if (any(is.na(simind))) {
warning("some simulated values were NA and were removed")
nsimul <- nsimul - rowSums(is.na(simind))
Modified: pkg/vegan/R/ordiareatest.R
===================================================================
--- pkg/vegan/R/ordiareatest.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/ordiareatest.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -16,6 +16,7 @@
function(ord, groups, area = c("hull", "ellipse"), permutations = 999,
parallel = getOption("mc.cores"), ...)
{
+ EPS <- sqrt(.Machine$double.eps)
## Function to find area
area <- match.arg(area)
areafun <- if (area == "hull") ordihull else ordiellipse
@@ -47,7 +48,7 @@
} else {
areas <- sapply(1:permutations, function(i, ...) pfun(perm[i,], ...))
}
- signif <- (rowSums(areas <= obs) + 1)/(nperm + 1)
+ signif <- (rowSums(areas <= obs + EPS) + 1)/(nperm + 1)
out <- list("areas" = obs, "pvalues" = signif, "permutations" = areas,
nperm = nperm, control = attr(perm, "control"), "kind" = area)
class(out) <- "ordiareatest"
Modified: pkg/vegan/R/permutest.betadisper.R
===================================================================
--- pkg/vegan/R/permutest.betadisper.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/permutest.betadisper.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
permutations = 999,
parallel = getOption("mc.cores"), ...)
{
+ EPS <- sqrt(.Machine$double.eps) # for P-value comparisons
t.statistic <- function(x, y) {
m <- length(x)
n <- length(y)
@@ -104,8 +105,8 @@
## Process results
F0 <- summary(mod)$fstatistic[1]
- Fstats <- round(Pstats[, 1], 12) # allow empty dim to be dropped
- statistic <- F0 <- round(F0, 12)
+ Fstats <- Pstats[, 1] # allow empty dim to be dropped
+ statistic <- F0
names(statistic) <- "Overall (F)"
## pairwise comparisons
@@ -113,13 +114,12 @@
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)
+ Tstats <- Pstats[, -1, drop = FALSE]
statistic <- c(statistic, T0)
}
## compute permutation p-value
- pval <- (sum(Fstats >= F0) + 1) / (length(Fstats) + 1)
+ pval <- (sum(Fstats >= F0 - EPS) + 1) / (length(Fstats) + 1)
if(pairwise) {
df <- apply(combin, 2, function(z) {
Modified: pkg/vegan/R/permutest.cca.R
===================================================================
--- pkg/vegan/R/permutest.cca.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/permutest.cca.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -139,10 +139,6 @@
num <- tmp[,1]
den <- tmp[,2]
F.perm <- tmp[,3]
- ## Round to avoid arbitrary ordering of statistics due to
- ## numerical inaccuracy
- F.0 <- round(F.0, 12)
- F.perm <- round(F.perm, 12)
Call <- match.call()
Call[[1]] <- as.name("permutest")
sol <- list(call = Call, testcall = x$call, model = model,
Modified: pkg/vegan/R/print.permutest.cca.R
===================================================================
--- pkg/vegan/R/print.permutest.cca.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/print.permutest.cca.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -1,10 +1,11 @@
`print.permutest.cca` <-
function (x, ...)
{
+ EPS <- sqrt(.Machine$double.eps)
cat("\nPermutation test for", x$method, "\n\n")
cat(howHead(x$control), "\n")
writeLines(strwrap(pasteCall(x$testcall)))
- Pval <- (sum(x$F.perm >= x$F.0) + 1)/(x$nperm + 1)
+ Pval <- (sum(x$F.perm >= x$F.0 - EPS) + 1)/(x$nperm + 1)
cat("Permutation test for ")
if (x$first)
cat("first constrained eigenvalue\n")
Modified: pkg/vegan/R/protest.R
===================================================================
--- pkg/vegan/R/protest.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/protest.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function (X, Y, scores = "sites", permutations = how(nperm = 999),
...)
{
+ EPS <- sqrt(.Machine$double.eps)
X <- scores(X, display = scores, ...)
Y <- scores(Y, display = scores, ...)
## Centre and normalize X & Y here so that the permutations will
@@ -34,7 +35,7 @@
perm <- sapply(seq_len(np),
function(i, ...) procr(X, Y[permutations[i,],]))
- Pval <- (sum(perm >= sol$t0) + 1)/(np + 1)
+ Pval <- (sum(perm >= sol$t0 - EPS) + 1)/(np + 1)
sol$t <- perm
sol$signif <- Pval
Modified: pkg/vegan/R/simper.R
===================================================================
--- pkg/vegan/R/simper.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/simper.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -2,6 +2,7 @@
function(comm, group, permutations = 0, trace = FALSE,
parallel = getOption("mc.cores"), ...)
{
+ EPS <- sqrt(.Machine$double.eps)
if (any(rowSums(comm, na.rm = TRUE) == 0))
warning("you have empty rows: results may be meaningless")
pfun <- function(x, comm, comp, i, contrp) {
@@ -78,7 +79,7 @@
perm.contr <- sapply(1:nperm, function(d)
pfun(d, comm, comp, i, contrp))
}
- p <- (rowSums(apply(perm.contr, 2, function(x) x >= average)) + 1) / (nperm + 1)
+ p <- (rowSums(apply(perm.contr, 2, function(x) x >= average - EPS)) + 1) / (nperm + 1)
}
else {
p <- NULL
Modified: pkg/vegan/R/vectorfit.R
===================================================================
--- pkg/vegan/R/vectorfit.R 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/R/vectorfit.R 2015-08-28 08:22:19 UTC (rev 2957)
@@ -1,6 +1,7 @@
"vectorfit" <-
function (X, P, permutations = 0, strata = NULL, w, ...)
{
+ EPS <- sqrt(.Machine$double.eps)
if (missing(w) || is.null(w))
w <- 1
if (length(w) == 1)
@@ -47,7 +48,7 @@
## permutations are the matrix columns and variables are rows
if (!is.matrix(permstore))
permstore <- matrix(permstore, ncol=permutations)
- permstore <- sweep(permstore, 1, r, ">=")
+ permstore <- sweep(permstore, 1, r - EPS, ">=")
validn <- rowSums(is.finite(permstore))
pvals <- (rowSums(permstore, na.rm = TRUE) + 1)/(validn + 1)
}
Modified: pkg/vegan/inst/NEWS.Rd
===================================================================
--- pkg/vegan/inst/NEWS.Rd 2015-08-28 07:40:29 UTC (rev 2956)
+++ pkg/vegan/inst/NEWS.Rd 2015-08-28 08:22:19 UTC (rev 2957)
@@ -7,6 +7,20 @@
\subsection{BUG FIXES}{
\itemize{
+ \item Permutation tests did not always correctly recognize ties
+ with the observed statistic and this could result in too low
+ \eqn{P}-values. This would happen in particular when all predictor
+ variables were factors (classes). The changes concern functions
+ \code{adonis}, \code{anosim}, \code{anova} and \code{permutest}
+ functions for \code{cca}, \code{rda} and \code{capscale},
+ \code{permutest} for \code{betadisper}, \code{envfit},
+ \code{mantel} and \code{mantel.partial}, \code{mrpp}, \code{mso},
+ \code{oecosimu}, \code{ordiareatest}, \code{protest} and
+ \code{simper}. This also fixes issues
+ \href{https://github.com/vegandevs/vegan/issues/120}{#120} and
+ \href{https://github.com/vegandevs/vegan/issues/132}{#132} in
+ GitHub.
+
\item Automated model building in constrained ordination
(\code{cca}, \code{rda}, \code{capscale}) with \code{step},
\code{ordistep} and \code{ordiR2step} could fail if there were
More information about the Vegan-commits
mailing list