[Vegan-commits] r1218 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jun 6 11:13:01 CEST 2010
Author: jarioksa
Date: 2010-06-06 11:13:01 +0200 (Sun, 06 Jun 2010)
New Revision: 1218
Modified:
branches/1.17/R/RsquareAdj.R
branches/1.17/R/adonis.R
branches/1.17/R/cca.formula.R
branches/1.17/R/mantel.R
branches/1.17/R/mantel.correlog.R
branches/1.17/R/mantel.partial.R
branches/1.17/R/ordiresids.R
branches/1.17/R/permutest.cca.R
branches/1.17/R/print.permutest.cca.R
branches/1.17/R/rda.formula.R
branches/1.17/R/scores.cca.R
branches/1.17/R/scores.rda.R
branches/1.17/inst/ChangeLog
branches/1.17/man/mantel.correlog.Rd
Log:
merge fixes and major enhancements to 1.17-3
Modified: branches/1.17/R/RsquareAdj.R
===================================================================
--- branches/1.17/R/RsquareAdj.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/RsquareAdj.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -8,10 +8,10 @@
`RsquareAdj.default` <-
function(x, n, m, ...)
{
- if (m >= (n-1))
- NA
- else
- 1 - (1-x)*(n-1)/(n-m-1)
+ r2 <- 1 - (1-x)*(n-1)/(n-m-1)
+ if (any(na <- m >= n-1))
+ r2[na] <- NA
+ r2
}
## Use this with rda() results
@@ -19,7 +19,7 @@
function(x, ...)
{
R2 <- x$CCA$tot.chi/x$tot.chi
- m <- x$CCA$rank
+ m <- x$CCA$qrank
n <- nrow(x$CCA$u)
if (is.null(x$pCCA))
radj <- RsquareAdj(R2, n, m)
@@ -33,7 +33,7 @@
function(x, ...)
{
R2 <- x$CCA$tot.chi/x$tot.chi
- m <- x$CCA$rank
+ m <- x$CCA$qrank
n <- nrow(x$CCA$u)
radj <- NA
list(r.squared = R2, adj.r.squared = radj)
Modified: branches/1.17/R/adonis.R
===================================================================
--- branches/1.17/R/adonis.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/adonis.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -84,14 +84,17 @@
f.test(H.s[[i]], G[p[,j],p[,j]], I, df.Exp[i], df.Res, H.snterm)
} )
})
-
+ ## Round to avoid arbitrary P-values with tied data
+ f.perms <- round(f.perms, 12)
+ F.Mod <- round(F.Mod, 12)
SumsOfSqs = c(SS.Exp.each, SS.Res, sum(SS.Exp.each) + SS.Res)
tab <- data.frame(Df = c(df.Exp, df.Res, n-1),
SumsOfSqs = SumsOfSqs,
MeanSqs = c(SS.Exp.each/df.Exp, SS.Res/df.Res, NA),
F.Model = c(F.Mod, NA,NA),
R2 = SumsOfSqs/SumsOfSqs[length(SumsOfSqs)],
- P = c((rowSums(t(f.perms) > F.Mod)+1)/(permutations+1), NA, NA))
+ P = c((rowSums(t(f.perms) >= F.Mod)+1)/(permutations+1),
+ NA, NA))
rownames(tab) <- c(attr(attr(rhs.frame, "terms"), "term.labels")[u.grps],
"Residuals", "Total")
colnames(tab)[ncol(tab)] <- "Pr(>F)"
Modified: branches/1.17/R/cca.formula.R
===================================================================
--- branches/1.17/R/cca.formula.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/cca.formula.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -14,10 +14,13 @@
sol$rowsum)
if (!is.null(sol$CCA$alias))
sol$CCA$centroids <- unique(sol$CCA$centroids)
+ ## See that there really are centroids
if (!is.null(sol$CCA$centroids)) {
rs <- rowSums(sol$CCA$centroids^2)
sol$CCA$centroids <- sol$CCA$centroids[rs > 1e-04, ,
drop = FALSE]
+ if (length(sol$CCA$centroids) == 0)
+ sol$CCA$centroids <- NULL
}
sol$terms <- d$terms
sol$terminfo <- ordiTerminfo(d, d$modelframe)
Modified: branches/1.17/R/mantel.R
===================================================================
--- branches/1.17/R/mantel.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/mantel.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -10,10 +10,13 @@
if (permutations) {
N <- attributes(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
+ xmat <- as.matrix(xdis)
+ asdist <- row(xmat) > col(xmat)
for (i in 1:permutations) {
take <- permuted.index(N, strata)
- permvec <- as.vector(as.dist(as.matrix(xdis)[take,
- take]))
+ permvec <- (xmat[take, take])[asdist]
perm[i] <- cor(permvec, ydis, method = method)
}
signif <- (sum(perm >= statistic) + 1)/(permutations + 1)
Modified: branches/1.17/R/mantel.correlog.R
===================================================================
--- branches/1.17/R/mantel.correlog.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/mantel.correlog.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -1,13 +1,13 @@
-`mantel.correlog` <-
- function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
- cutoff=TRUE, r.type="pearson", nperm=999, mult="holm",
- progressive=TRUE)
-{
+`mantel.correlog` <-
+ function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
+ cutoff=TRUE, r.type="pearson", nperm=999, mult="holm",
+ progressive=TRUE)
+{
r.type <- match.arg(r.type, c("pearson", "spearman", "kendall"))
mult <- match.arg(mult, c("sidak", p.adjust.methods))
-
+
epsilon <- .Machine$double.eps
- D.eco = as.matrix(D.eco)
+ D.eco <- as.matrix(D.eco)
## Geographic distance matrix
if(!is.null(D.geo)) {
@@ -16,7 +16,7 @@
D.geo <- as.matrix(D.geo)
} else {
if(is.null(XY)) {
- stop("You did not provided a geographic distance matrix nor a list of site coordinates.")
+ stop("You did not provide a geographic distance matrix nor a list of site coordinates")
} else {
D.geo <- as.matrix(dist(XY))
}
@@ -24,46 +24,49 @@
n <- nrow(D.geo)
if(n != nrow(D.eco))
- stop("Numbers of objects in D.eco and D.geo are not equal.")
+ stop("Numbers of objects in D.eco and D.geo are not equal")
n.dist <- n*(n-1)/2
vec.D <- as.vector(as.dist(D.geo))
vec.DD <- as.vector(D.geo)
-
- ## Number of classes
- if(n.class > 0) {
- if(!is.null(break.pts))
+
+ ## Number of classes and breakpoints
+
+ if(!is.null(break.pts)) {
+ ## Use the list of break points
+ if(n.class > 0)
stop("You provided both a number of classes and a list of break points. Which one should the function use?")
+ n.class = length(break.pts) - 1
+
} else {
- if(is.null(break.pts)) {
- ## Use Sturge's rule to determine the number of classes
- n.class <- round(1 + log(n.dist, base=2))
- start.pt <- min(vec.D)
- end.pt <- max(vec.D)
- width <- (end.pt - start.pt)/n.class
- break.pts <- vector(length=n.class+1)
- break.pts[n.class+1] <- end.pt
- for(i in 1:n.class) {
- break.pts[i] <- start.pt + width*(i-1) }
- } else {
- ## Use the list of break points
- n.class <- length(break.pts) - 1
+ ## No breakpoints have been provided: equal-width classes
+ if(n.class == 0) {
+ ## Use Sturges rule to determine the number of classes
+ n.class <- ceiling(1 + log(n.dist, base=2))
}
+ ## Compute the breakpoints from n.class
+ start.pt <- min(vec.D)
+ end.pt <- max(vec.D)
+ width <- (end.pt - start.pt)/n.class
+ break.pts <- vector(length=n.class+1)
+ break.pts[n.class+1] <- end.pt
+ for(i in 1:n.class)
+ break.pts[i] <- start.pt + width*(i-1)
}
+
half.cl <- n.class %/% 2
-
+
## Move the first breakpoint a little bit to the left
- break.pts[1] = break.pts[1] - epsilon
-
+ break.pts[1] <- break.pts[1] - epsilon
+
## Find the break points and the class indices
class.ind <- break.pts[1:n.class] +
(0.5*(break.pts[2:(n.class+1)]-break.pts[1:n.class]))
-
+
## Create the matrix of distance classes
vec2 <- vector(length=n^2)
- for(i in 1:n^2) {
- vec2[i] <- min( which(break.pts >= vec.DD[i]) ) - 1
- }
-
+ for(i in 1:n^2)
+ vec2[i] <- min( which(break.pts >= vec.DD[i]) ) - 1
+
## Start assembling the vectors of results
class.index <- NA
n.dist <- NA
@@ -107,7 +110,7 @@
mantel.res <- cbind(class.index, n.dist, mantel.r, mantel.p)
mantel.res <- mantel.res[-1,]
-
+
## Note: vector 'mantel.p' starts with a NA value
mantel.p <- mantel.p[-1]
n.tests <- length(which(!is.na(mantel.p)))
@@ -120,8 +123,8 @@
if(progressive) {
p.corr <- mantel.p[1]
if(mult == "sidak") {
- for(j in 2:n.tests) {
- p.corr <- c(p.corr, 1-(1-mantel.p[j])^j) }
+ for(j in 2:n.tests)
+ p.corr <- c(p.corr, 1-(1-mantel.p[j])^j)
} else {
for(j in 2:n.tests) {
temp <- p.adjust(mantel.p[1:j], method=mult)
@@ -150,4 +153,3 @@
class(res) <- "mantel.correlog"
res
}
-
Modified: branches/1.17/R/mantel.partial.R
===================================================================
--- branches/1.17/R/mantel.partial.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/mantel.partial.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -17,10 +17,11 @@
if (permutations) {
N <- attributes(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)
- permvec <- as.vector(as.dist(as.matrix(xdis)[take,
- 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)
Modified: branches/1.17/R/ordiresids.R
===================================================================
--- branches/1.17/R/ordiresids.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/ordiresids.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -8,6 +8,12 @@
stop("function is only available for constrained ordination")
fit <- fitted(x, type = residuals)
res <- residuals(x, type = residuals)
+ ## remove the effects of row weights in CA
+ if (!inherits(x, "rda")) {
+ sqr <- sqrt(x$rowsum)
+ fit <- sweep(fit, 1, sqr, "*")
+ res <- sweep(res, 1, sqr, "*")
+ }
colnam <- rep(colnames(fit), each=nrow(fit))
rownam <- rep(rownames(fit), ncol(fit))
df <- data.frame(Fitted = as.vector(fit), Residuals = as.vector(res))
Modified: branches/1.17/R/permutest.cca.R
===================================================================
--- branches/1.17/R/permutest.cca.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/permutest.cca.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -51,26 +51,36 @@
else E <- x$CA$Xbar
if (isPartial && model == "direct")
E <- E + Y.Z
+ ## Save dimensions
N <- nrow(E)
+ if (isCCA) {
+ Xcol <- ncol(X)
+ if (isPartial)
+ Zcol <- ncol(Z)
+ }
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
for (i in 1:permutations) {
take <- permuted.index(N, strata)
Y <- E[take, ]
+ if (isCCA)
+ wtake <- w[take]
if (isPartial) {
if (isCCA) {
- wm <- colSums(sweep(Z, 1, w[take], "*"))
- XZ <- sweep(Z, 2, wm, "-")
- XZ <- sweep(XZ, 1, sqrt(w[take]), "*")
+ XZ <- .C("wcentre", x = as.double(Z), as.double(wtake),
+ as.integer(N), as.integer(Zcol),
+ PACKAGE = "vegan")$x
+ dim(XZ) <- c(N, Zcol)
QZ <- qr(XZ)
}
Y <- qr.resid(QZ, Y)
}
if (isCCA) {
- wm <- colSums(sweep(X, 1, w[take], "*"))
- XY <- sweep(X, 2, wm, "-")
- XY <- sweep(XY, 1, sqrt(w[take]), "*")
+ XY <- .C("wcentre", x = as.double(X), as.double(wtake),
+ as.integer(N), as.integer(Xcol),
+ PACKAGE = "vegan")$x
+ dim(XY) <- c(N, Xcol)
Q <- qr(XY)
}
tmp <- qr.fitted(Q, Y)
Modified: branches/1.17/R/print.permutest.cca.R
===================================================================
--- branches/1.17/R/print.permutest.cca.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/print.permutest.cca.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -3,7 +3,7 @@
{
cat("\nPermutation test for", x$method, "\n\n")
writeLines(strwrap(pasteCall(x$call)))
- Pval <- sum(x$F.perm >= x$F.0)/x$nperm
+ Pval <- (sum(x$F.perm >= x$F.0) + 1)/(x$nperm + 1)
cat("Permutation test for ")
if (x$first)
cat("first constrained eigenvalue\n")
Modified: branches/1.17/R/rda.formula.R
===================================================================
--- branches/1.17/R/rda.formula.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/rda.formula.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -18,6 +18,8 @@
rs <- rowSums(sol$CCA$centroids^2)
sol$CCA$centroids <- sol$CCA$centroids[rs > 1e-04, ,
drop = FALSE]
+ if (length(sol$CCA$centroids) == 0)
+ sol$CCA$centroids <- NULL
}
sol$terms <- d$terms
sol$terminfo <- ordiTerminfo(d, d$modelframe)
Modified: branches/1.17/R/scores.cca.R
===================================================================
--- branches/1.17/R/scores.cca.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/scores.cca.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -90,6 +90,14 @@
sol$centroids <- cn
}
}
+ ## Take care that scores have names
+ for (i in 1:length(sol)) {
+ if (is.matrix(sol[[i]]))
+ rownames(sol[[i]]) <-
+ rownames(sol[[i]], do.NULL = FALSE,
+ prefix = substr(names(sol)[i], 1, 3))
+ }
+ ## Only one type of scores: return a matrix instead of a list
if (length(sol) == 1)
sol <- sol[[1]]
return(sol)
Modified: branches/1.17/R/scores.rda.R
===================================================================
--- branches/1.17/R/scores.rda.R 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/R/scores.rda.R 2010-06-06 09:13:01 UTC (rev 1218)
@@ -82,7 +82,15 @@
}
sol$centroids <- cn
}
- }
+ }
+ ## Take care that scores have names
+ for (i in 1:length(sol)) {
+ if (is.matrix(sol[[i]]))
+ rownames(sol[[i]]) <-
+ rownames(sol[[i]], do.NULL = FALSE,
+ prefix = substr(names(sol)[i], 1, 3))
+ }
+ ## Only one type of scores: return a matrix instead of a list
if (length(sol) == 1)
sol <- sol[[1]]
attr(sol, "const") <- const
Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/inst/ChangeLog 2010-06-06 09:13:01 UTC (rev 1218)
@@ -4,6 +4,28 @@
Version 1.17-3 (not yet released)
+ * merged r1210:1213: mantel and mantel.partial speed-ups.
+
+ * merge r1200: cca, rda failed when Condition() was a factor, but
+ constraints had no factors.
+
+ * merge r1190: RsquareAdj.default handles vector arguments.
+
+ * merge r1188: ordiresids de-weights *CA.
+
+ * merge r1182 (r1184 for Rd) r1187: mantel.correlog upgrades.
+
+ * merge r1181: RsquareAdj.rda used wrong df in rank-deficit
+ models.
+
+ * merge r1179: tie handling in adonis.
+
+ * merge r1177: print.permutest.cca uses correct order statistic.
+
+ * merge r1176: permutest.cca speed-up in CCA re-weighting.
+
+ * merge r1174: scores.cca/rda always have names.
+
* betadisper: 'type = "median"', the default, was not computing
the spatial median on the real and imaginary axes separately.
Reported by Marek Omelka.
Modified: branches/1.17/man/mantel.correlog.Rd
===================================================================
--- branches/1.17/man/mantel.correlog.Rd 2010-06-03 18:29:56 UTC (rev 1217)
+++ branches/1.17/man/mantel.correlog.Rd 2010-06-06 09:13:01 UTC (rev 1218)
@@ -28,11 +28,12 @@
\item{XY}{ A file of Cartesian geographic coordinates of the
points. Default: \code{XY=NULL}. }
- \item{n.class}{ Number of classes. If \code{n.class=0}, the Sturge
+ \item{n.class}{ Number of classes. If \code{n.class=0}, the Sturges
equation will be used unless break points are provided. }
\item{break.pts}{ Vector containing the break points of the distance
- distribution. Default: \code{break.pts=NULL}. }
+ distribution. Provide (n.class+1) breakpoints, that is, a list with
+ a beginning and an ending point. Default: \code{break.pts=NULL}. }
\item{cutoff}{ For the second half of the distance classes,
\code{cutoff = TRUE} limits the correlogram to the distance classes
@@ -102,16 +103,18 @@
the program. }
\item{mult }{The name of the correction for multiple testing. No
- correction: \code{mult="none"}. } #
+ correction: \code{mult="none"}. }
\item{progressive }{A logical (\code{TRUE}, \code{FALSE}) value
indicating whether or not a progressive correction for multiple
- testing was requested. } \item{n.tests }{The number of distance
- classes for which Mantel tests have been computed and tested for
- significance. }
+ testing was requested. }
-\item{call }{The function call. } }
+ \item{n.tests }{The number of distance classes for which Mantel
+ tests have been computed and tested for significance. }
+ \item{call }{The function call. }
+}
+
\author{ Pierre Legendre, Universite de Montreal }
\references{
@@ -128,32 +131,38 @@
Sokal, R. R. 1986. Spatial data analysis and historical
processes. 29-43 in: E. Diday et al. [eds.] Data analysis and
- informatics, IV. North-Holland, Amsterdam. }
+ informatics, IV. North-Holland, Amsterdam.
+
+ Sturges, H. A. 1926. The choice of a class interval. Journal of the
+ American Statistical Association 21: 65–66. }
\examples{
-# Mite data from "vegan"
+# Mite data available in "vegan"
data(mite)
data(mite.xy)
mite.hel <- decostand(mite, "hellinger")
-mite.hel.D <- dist(mite.hel)
-mite.correlog <- mantel.correlog(mite.hel.D, XY=mite.xy, nperm=99)
+# Detrend the species data by regression on the site coordinates
+mite.hel.resid <- resid(lm(as.matrix(mite.hel) ~ ., data=mite.xy))
+
+# Compute the detrended species distance matrix
+mite.hel.D = dist(mite.hel.resid)
+
+# Compute Mantel correlogram with cutoff, Pearson statistic
+mite.correlog = mantel.correlog(mite.hel.D, XY=mite.xy, nperm=99)
summary(mite.correlog)
mite.correlog
+# or: print(mite.correlog)
+# or: print.mantel.correlog(mite.correlog)
plot(mite.correlog)
-mite.correlog2 <- mantel.correlog(mite.hel.D, XY=mite.xy, cutoff=FALSE,
+# Compute Mantel correlogram without cutoff, Spearman statistic
+mite.correlog2 = mantel.correlog(mite.hel.D, XY=mite.xy, cutoff=FALSE,
r.type="spearman", nperm=99)
summary(mite.correlog2)
mite.correlog2
plot(mite.correlog2)
-## Mite correlogram after spatially detrending the mite data
-mite.h.det <- resid(lm(as.matrix(mite.hel.D) ~ ., data=mite.xy))
-mite.correlog3 <- mantel.correlog(mite.h.det, XY=mite.xy, nperm=99)
-mite.correlog3
-plot(mite.correlog3)
-
}
\keyword{ multivariate }
\ No newline at end of file
More information about the Vegan-commits
mailing list