[Vegan-commits] r847 - in branches/1.15: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 1 10:36:37 CEST 2009
Author: jarioksa
Date: 2009-06-01 10:36:37 +0200 (Mon, 01 Jun 2009)
New Revision: 847
Modified:
branches/1.15/R/adonis.R
branches/1.15/R/anosim.R
branches/1.15/R/factorfit.R
branches/1.15/R/mantel.R
branches/1.15/R/mantel.partial.R
branches/1.15/R/mrpp.R
branches/1.15/R/plot.anosim.R
branches/1.15/R/print.anosim.R
branches/1.15/R/print.factorfit.R
branches/1.15/R/print.mantel.R
branches/1.15/R/print.mrpp.R
branches/1.15/R/print.permutest.cca.R
branches/1.15/R/print.protest.R
branches/1.15/R/print.vectorfit.R
branches/1.15/R/vectorfit.R
branches/1.15/inst/ChangeLog
branches/1.15/man/anosim.Rd
branches/1.15/man/envfit.Rd
branches/1.15/man/mantel.Rd
branches/1.15/man/mrpp.Rd
Log:
put observed statistic among permutations: merged rev 837:842 and 845 to branches/1.15
Modified: branches/1.15/R/adonis.R
===================================================================
--- branches/1.15/R/adonis.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/adonis.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -8,11 +8,6 @@
## variables. data is the data frame from which A, B, and C would
## be drawn.
TOL <- 1e-7
- ## Set no. of permutations to x-1 if x is an even hundred
- if (permutations %% 100 == 0) {
- permutations <- permutations - 1
- warning("Setting no. of permutations to ", permutations)
- }
Terms <- terms(formula, data = data)
lhs <- formula[[2]]
lhs <- eval(lhs, data, parent.frame()) # to force evaluation
Modified: branches/1.15/R/anosim.R
===================================================================
--- branches/1.15/R/anosim.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/anosim.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,5 +1,5 @@
"anosim" <-
- function (dis, grouping, permutations = 1000, strata)
+ function (dis, grouping, permutations = 999, strata)
{
x <- as.dist(dis)
sol <- c(call = match.call())
@@ -28,7 +28,7 @@
tmp.ave <- tapply(x.rank, tmp.within, mean)
perm[i] <- -diff(tmp.ave)/div
}
- p.val <- sum(perm >= statistic)/permutations
+ p.val <- (1 + sum(perm >= statistic))/(1 + permutations)
sol$signif <- p.val
sol$perm <- perm
}
Modified: branches/1.15/R/factorfit.R
===================================================================
--- branches/1.15/R/factorfit.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/factorfit.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -45,7 +45,7 @@
var = double(1), PACKAGE = "vegan")$var
tmp[i] <- 1 - invar/totvar
}
- pval.this <- sum(tmp > r.this)/permutations
+ pval.this <- (sum(tmp > r.this) + 1)/(permutations + 1)
pval <- c(pval, pval.this)
}
}
Modified: branches/1.15/R/mantel.R
===================================================================
--- branches/1.15/R/mantel.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/mantel.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,33 +1,33 @@
"mantel" <-
- function (xdis, ydis, method = "pearson", permutations = 1000,
+ function (xdis, ydis, method = "pearson", permutations = 999,
strata)
{
- xdis <- as.dist(xdis)
- ydis <- as.vector(as.dist(ydis))
- tmp <- cor.test(as.vector(xdis), ydis, method = method)
- statistic <- as.numeric(tmp$estimate)
- variant <- tmp$method
- if (permutations) {
- N <- attributes(xdis)$Size
- perm <- rep(0, permutations)
- for (i in 1:permutations) {
- take <- permuted.index(N, strata)
- permvec <- as.vector(as.dist(as.matrix(xdis)[take,
- take]))
- perm[i] <- cor(permvec, ydis, method = method)
+ xdis <- as.dist(xdis)
+ ydis <- as.vector(as.dist(ydis))
+ tmp <- cor.test(as.vector(xdis), ydis, method = method)
+ statistic <- as.numeric(tmp$estimate)
+ variant <- tmp$method
+ if (permutations) {
+ N <- attributes(xdis)$Size
+ perm <- rep(0, permutations)
+ for (i in 1:permutations) {
+ take <- permuted.index(N, strata)
+ permvec <- as.vector(as.dist(as.matrix(xdis)[take,
+ take]))
+ perm[i] <- cor(permvec, ydis, method = method)
+ }
+ signif <- (sum(perm >= statistic) + 1)/(permutations + 1)
+ }
+ else {
+ signif <- NA
+ perm <- NULL
}
- signif <- sum(perm >= statistic)/permutations
- }
- else {
- signif <- NA
- perm <- NULL
- }
- res <- list(call = match.call(), method = variant, statistic = statistic,
- signif = signif, perm = perm, permutations = permutations)
- if (!missing(strata)) {
- res$strata <- deparse(substitute(strata))
- res$stratum.values <- strata
- }
- class(res) <- "mantel"
- res
+ res <- list(call = match.call(), method = variant, statistic = statistic,
+ signif = signif, perm = perm, permutations = permutations)
+ if (!missing(strata)) {
+ res$strata <- deparse(substitute(strata))
+ res$stratum.values <- strata
+ }
+ class(res) <- "mantel"
+ res
}
Modified: branches/1.15/R/mantel.partial.R
===================================================================
--- branches/1.15/R/mantel.partial.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/mantel.partial.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,43 +1,43 @@
"mantel.partial" <-
- function (xdis, ydis, zdis, method = "pearson", permutations = 1000,
+ function (xdis, ydis, zdis, method = "pearson", permutations = 999,
strata)
{
- part.cor <- function(rxy, rxz, ryz) {
- (rxy - rxz * ryz)/sqrt(1-rxz*rxz)/sqrt(1-ryz*ryz)
- }
- xdis <- as.dist(xdis)
- ydis <- as.vector(as.dist(ydis))
- zdis <- as.vector(as.dist(zdis))
- rxy <- cor.test(as.vector(xdis), ydis, method = method)
- rxz <- cor(as.vector(xdis), zdis, method = method)
- ryz <- cor(ydis, zdis, method = method)
- variant <- rxy$method
- rxy <- rxy$estimate
- statistic <- part.cor(rxy, rxz, ryz)
- if (permutations) {
- N <- attributes(xdis)$Size
- perm <- rep(0, permutations)
- for (i in 1:permutations) {
- take <- permuted.index(N, strata)
- permvec <- as.vector(as.dist(as.matrix(xdis)[take,
- take]))
- rxy <- cor(permvec, ydis, method = method)
- rxz <- cor(permvec, zdis, method = method)
- perm[i] <- part.cor(rxy, rxz, ryz)
+ part.cor <- function(rxy, rxz, ryz) {
+ (rxy - rxz * ryz)/sqrt(1-rxz*rxz)/sqrt(1-ryz*ryz)
}
- signif <- sum(perm >= statistic)/permutations
- }
- else {
- signif <- NA
- perm <- NULL
- }
- res <- list(call = match.call(), method = variant, statistic = statistic,
- signif = signif, perm = perm, permutations = permutations)
- if (!missing(strata)) {
- res$strata <- deparse(substitute(strata))
- res$stratum.values <- strata
- }
- class(res) <- c("mantel.partial", "mantel")
- res
+ xdis <- as.dist(xdis)
+ ydis <- as.vector(as.dist(ydis))
+ zdis <- as.vector(as.dist(zdis))
+ rxy <- cor.test(as.vector(xdis), ydis, method = method)
+ rxz <- cor(as.vector(xdis), zdis, method = method)
+ ryz <- cor(ydis, zdis, method = method)
+ variant <- rxy$method
+ rxy <- rxy$estimate
+ statistic <- part.cor(rxy, rxz, ryz)
+ if (permutations) {
+ N <- attributes(xdis)$Size
+ perm <- rep(0, permutations)
+ for (i in 1:permutations) {
+ take <- permuted.index(N, strata)
+ permvec <- as.vector(as.dist(as.matrix(xdis)[take,
+ take]))
+ rxy <- cor(permvec, ydis, method = method)
+ rxz <- cor(permvec, zdis, method = method)
+ perm[i] <- part.cor(rxy, rxz, ryz)
+ }
+ signif <- (sum(perm >= statistic)+1)/(permutations + 1)
+ }
+ else {
+ signif <- NA
+ perm <- NULL
+ }
+ res <- list(call = match.call(), method = variant, statistic = statistic,
+ signif = signif, perm = perm, permutations = permutations)
+ if (!missing(strata)) {
+ res$strata <- deparse(substitute(strata))
+ res$stratum.values <- strata
+ }
+ class(res) <- c("mantel.partial", "mantel")
+ res
}
Modified: branches/1.15/R/mrpp.R
===================================================================
--- branches/1.15/R/mrpp.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/mrpp.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,5 +1,5 @@
"mrpp" <-
-function (dat, grouping, permutations = 1000, distance = "euclidean",
+function (dat, grouping, permutations = 999, distance = "euclidean",
weight.type = 1, strata)
{
classmean <- function(ind, dmat, indls) {
Modified: branches/1.15/R/plot.anosim.R
===================================================================
--- branches/1.15/R/plot.anosim.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/plot.anosim.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -5,7 +5,7 @@
...)
title(title)
if (x$permutations) {
- pval <- format.pval(x$signif, eps=1/x$permutations)
+ pval <- format.pval(x$signif)
} else {
pval <- "not assessed"
}
Modified: branches/1.15/R/print.anosim.R
===================================================================
--- branches/1.15/R/print.anosim.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.anosim.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -8,7 +8,7 @@
cat(formatC(x$statistic, digits = digits), "\n")
nperm <- x$permutations
if (nperm) {
- cat(" Significance:", format.pval(x$signif, eps = 1/nperm),
+ cat(" Significance:", format.pval(x$signif),
"\n\n")
cat("Based on ", nperm, " permutations")
}
Modified: branches/1.15/R/print.factorfit.R
===================================================================
--- branches/1.15/R/print.factorfit.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.factorfit.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,15 +1,12 @@
"print.factorfit" <-
function (x, ...)
{
- if (x$permutations)
- eps <- 1/x$permutations
- else eps <- .Machine$double.eps
cat("Centroids:\n")
printCoefmat(x$centroids, tst.ind = 1:ncol(x$centroids), na.print = "", ...)
cat("\nGoodness of fit:\n")
out <- cbind(r2 = x$r, "Pr(>r)" = x$pvals)
if (x$permutations) {
- printCoefmat(out, has.Pvalue = TRUE, eps.Pvalue = eps, ...)
+ printCoefmat(out, has.Pvalue = TRUE, ...)
cat("P values based on", x$permutations, "permutations")
if (!is.null(x$strata))
cat(", stratified within", x$strata)
Modified: branches/1.15/R/print.mantel.R
===================================================================
--- branches/1.15/R/print.mantel.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.mantel.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -11,7 +11,7 @@
cat(formatC(x$statistic, digits = digits), "\n")
nperm <- x$permutations
if (nperm) {
- cat(" Significance:", format.pval(x$signif, eps = 1/nperm),
+ cat(" Significance:", format.pval(x$signif),
"\n\n")
out <- quantile(x$perm, c(0.9, 0.95, 0.975, 0.99))
cat("Empirical upper confidence limits of r:\n")
Modified: branches/1.15/R/print.mrpp.R
===================================================================
--- branches/1.15/R/print.mrpp.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.mrpp.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -22,7 +22,7 @@
formatC(x$E.delta),"\n\n")
nperm <- x$permutations
if (nperm) {
- cat("Significance of delta:", format.pval(x$Pvalue, eps = 1/nperm),
+ cat("Significance of delta:", format.pval(x$Pvalue),
"\n")
cat("Based on ", nperm, " permutations")
}
Modified: branches/1.15/R/print.permutest.cca.R
===================================================================
--- branches/1.15/R/print.permutest.cca.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.permutest.cca.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -10,7 +10,7 @@
else
cat("all constrained eigenvalues\n")
cat("Pseudo-F:\t", x$F.0, "\n")
- cat("Significance:\t", format.pval(Pval, eps = 1/x$nperm),
+ cat("Significance:\t", format.pval(Pval),
"\n")
cat("Based on", x$nperm, "permutations under", x$model, "model")
if (!is.null(x$strata))
Modified: branches/1.15/R/print.protest.R
===================================================================
--- branches/1.15/R/print.protest.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.protest.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -6,7 +6,7 @@
cat("Correlation in a symmetric Procrustes rotation: ")
cat(formatC(x$t0, digits = digits), "\n")
cat("Significance: ")
- cat(format.pval(x$signif, eps = 1/x$permutations),"\n")
+ cat(format.pval(x$signif),"\n")
cat("Based on", x$permutations, "permutations")
if (!is.null(x$strata))
cat(", stratified within", x$strata)
Modified: branches/1.15/R/print.vectorfit.R
===================================================================
--- branches/1.15/R/print.vectorfit.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/print.vectorfit.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -1,11 +1,8 @@
"print.vectorfit" <-
function (x, ...)
{
- if (x$permutations)
- eps <- 1/x$permutations
- else eps <- .Machine$double.eps
out <- cbind(x$arrows, r2 = x$r, "Pr(>r)" = x$pvals)
- printCoefmat(out, eps.Pvalue = eps, na.print = "", ...)
+ printCoefmat(out, na.print = "", ...)
if (x$permutations) {
cat("P values based on", x$permutations, "permutations")
if (!is.null(x$strata))
Modified: branches/1.15/R/vectorfit.R
===================================================================
--- branches/1.15/R/vectorfit.R 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/R/vectorfit.R 2009-06-01 08:36:37 UTC (rev 847)
@@ -42,7 +42,7 @@
permstore[i, ] <- diag(cor(Hperm, take))^2
}
permstore <- sweep(permstore, 2, r, ">")
- pvals <- apply(permstore, 2, sum)/permutations
+ pvals <- (apply(permstore, 2, sum) + 1)/(permutations + 1)
}
else pvals <- NULL
sol <- list(arrows = heads, r = r, permutations = permutations,
Modified: branches/1.15/inst/ChangeLog
===================================================================
--- branches/1.15/inst/ChangeLog 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/inst/ChangeLog 2009-06-01 08:36:37 UTC (rev 847)
@@ -5,6 +5,12 @@
Version 1.15-3 (opened 28 May, 2009)
+ * merged r837:r842 and r846 that put the observed test statistic
+ among permutations like should be done with order statistics.
+ Concerns functions anosim, mantel, mantel.partial, mrpp, protest,
+ factorfit, vectorfit. General cleanup also concerns adonis and
+ permutest.cca.
+
* merged r832: wrong reference in metaMDS.Rd.
* merged r827: kendall.global: could get wrong counts of ties in
Modified: branches/1.15/man/anosim.Rd
===================================================================
--- branches/1.15/man/anosim.Rd 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/man/anosim.Rd 2009-06-01 08:36:37 UTC (rev 847)
@@ -11,7 +11,7 @@
of sampling units.
}
\usage{
-anosim(dis, grouping, permutations=1000, strata)
+anosim(dis, grouping, permutations=999, strata)
}
\arguments{
Modified: branches/1.15/man/envfit.Rd
===================================================================
--- branches/1.15/man/envfit.Rd 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/man/envfit.Rd 2009-06-01 08:36:37 UTC (rev 847)
@@ -160,7 +160,7 @@
data(varechem)
library(MASS)
ord <- metaMDS(varespec)
-(fit <- envfit(ord, varechem, perm = 1000))
+(fit <- envfit(ord, varechem, perm = 999))
scores(fit, "vectors")
plot(ord)
plot(fit)
@@ -169,7 +169,7 @@
## that arrows are scaled similarly in cca and envfit plots
ord <- cca(varespec ~ Al + P + K, varechem)
plot(ord, type="p")
-fit <- envfit(ord, varechem, perm = 1000, display = "lc")
+fit <- envfit(ord, varechem, perm = 999, display = "lc")
plot(fit, p.max = 0.05, col = "red")
## Class variables, formula interface, and displaying the
## inter-class variability with `ordispider'
Modified: branches/1.15/man/mantel.Rd
===================================================================
--- branches/1.15/man/mantel.Rd 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/man/mantel.Rd 2009-06-01 08:36:37 UTC (rev 847)
@@ -14,8 +14,8 @@
}
\usage{
-mantel(xdis, ydis, method="pearson", permutations=1000, strata)
-mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 1000,
+mantel(xdis, ydis, method="pearson", permutations=999, strata)
+mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 999,
strata)
}
Modified: branches/1.15/man/mrpp.Rd
===================================================================
--- branches/1.15/man/mrpp.Rd 2009-06-01 08:30:55 UTC (rev 846)
+++ branches/1.15/man/mrpp.Rd 2009-06-01 08:36:37 UTC (rev 847)
@@ -10,7 +10,7 @@
groups of sampling units. }
\usage{
-mrpp(dat, grouping, permutations = 1000, distance = "euclidean",
+mrpp(dat, grouping, permutations = 999, distance = "euclidean",
weight.type = 1, strata)
}
More information about the Vegan-commits
mailing list