[Vegan-commits] r339 - in pkg: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 7 15:54:12 CEST 2008
Author: jarioksa
Date: 2008-05-07 15:54:12 +0200 (Wed, 07 May 2008)
New Revision: 339
Added:
pkg/R/mso.R
pkg/R/plot.mso.R
pkg/man/mso.Rd
Modified:
pkg/inst/ChangeLog
Log:
added Helene Wagner's mso function for spatial partioning of rda and cca results
Added: pkg/R/mso.R
===================================================================
--- pkg/R/mso.R (rev 0)
+++ pkg/R/mso.R 2008-05-07 13:54:12 UTC (rev 339)
@@ -0,0 +1,76 @@
+mso <- function (object.cca, object.xy, grain = 1, round.up = FALSE,
+ permutations = FALSE)
+{
+ object <- object.cca
+ xy <- object.xy
+ N <- nrow(object$CA$Xbar)
+ if (class(object) == "rda")
+ N <- 1
+ require(vegan)
+ Dist <- dist(xy)
+ object$grain <- grain
+ if (round.up)
+ H <- ceiling(Dist/grain) * grain
+ else H <- round(Dist/grain) * grain
+ object$H <- H
+ H <- as.vector(H)
+ Dist <- sapply(split(Dist, H), mean)
+ object$vario <- data.frame(H = names(table(H)), Dist = Dist,
+ n = as.numeric(table(H)))
+ test <- list()
+ if (is.numeric(object$CCA$rank)) {
+ if (is.numeric(object$pCCA$rank)) {
+ test$pcca <- sapply(split(dist(object$pCCA$Fit)^2 *
+ N/2, H), mean)
+ test$cca <- sapply(split(dist(object$CCA$Xbar - object$CA$Xbar)^2 *
+ N/2, H), mean)
+ test$ca <- sapply(split(dist(object$CA$Xbar)^2 *
+ N/2, H), mean)
+ test$all.cond <- sapply(split(dist(object$CCA$Xbar)^2 *
+ N/2, H), mean)
+ test$se <- sqrt(sapply(split(dist(object$CCA$Xbar)^2 *
+ N/2, H), var)/object$vario$n)
+ object$vario <- cbind(object$vario, All = test$all.cond,
+ Sum = test$ca + test$cca, CA = test$ca,
+ CCA = test$cca, pCCA = test$pcca,
+ se = test$se)
+ } else {
+ test$all <- sapply(split(dist(object$CCA$Xbar)^2 *
+ N/2, H), mean)
+ test$cca <- sapply(split(dist(object$CCA$Xbar - object$CA$Xbar)^2 *
+ N/2, H), mean)
+ test$ca <- sapply(split(dist(object$CA$Xbar)^2 *
+ N/2, H), mean)
+ test$se <- sqrt(sapply(split(dist(object$CCA$Xbar)^2 *
+ N/2, H), var)/object$vario$n)
+ object$vario <- cbind(object$vario, All = test$all,
+ Sum = test$ca + test$cca, CA = test$ca,
+ CCA = test$cca, se = test$se)
+ }
+ } else {
+ test$ca <- sapply(split(dist(object$CA$Xbar)^2 * N/2,
+ H), mean)
+ object$vario <- cbind(object$vario, All = test$ca, CA = test$ca)
+ }
+ if (permutations) {
+ ##require(base)
+ object$H.test <- matrix(0, length(object$H), nrow(object$vario))
+ for (i in 1:nrow(object$vario)) {
+ object$H.test[, i] <- as.numeric(object$H == object$vario$H[i])
+ }
+ xdis <- dist(object$CA$Xbar)^2
+ N <- attributes(xdis)$Size
+ statistic <- abs(cor(as.vector(xdis), object$H.test))
+ perm <- matrix(0, length(statistic), permutations)
+ for (i in 1:permutations) {
+ take <- sample(N, N)
+ permvec <- as.vector(as.dist(as.matrix(xdis)[take,
+ take]))
+ perm[, i] <- abs(cor(permvec, object$H.test))
+ }
+ object$vario$CA.signif <- apply((perm >= matrix(statistic,
+ nrow(perm), ncol(perm)))/permutations, 1, sum)
+ }
+ class(object) <- c("mso", class(object))
+ object
+}
Added: pkg/R/plot.mso.R
===================================================================
--- pkg/R/plot.mso.R (rev 0)
+++ pkg/R/plot.mso.R 2008-05-07 13:54:12 UTC (rev 339)
@@ -0,0 +1,91 @@
+plot.mso <- function (x, alpha = 0.05, explained = FALSE, ...)
+{
+ object.cca <- x
+ if (is.data.frame(object.cca$vario)) {
+ object <- object.cca
+ vario <- object$vario
+ grain <- object$grain
+ z <- qnorm(alpha/2)
+ if (is.numeric(vario$CA.signif)) {
+ vario <- vario[, -ncol(vario)]
+ }
+ ymax <- max(vario[, -1:-3], na.rm = TRUE)
+ b <- ncol(vario) - 3
+ label <- c("", "", "", "Total variance", "Explained plus residual",
+ "Residual variance", "Explained variance", "Conditioned variance")
+ par(omi = c(0.5, 0.5, 0, 0))
+ if (is.numeric(object$CCA$rank)) {
+ if (!explained)
+ b <- b - 1
+ if (is.numeric(object$vario$se))
+ b <- b - 1
+ plot(vario$Dist, vario$All, type = "n", lty = 1,
+ pch = 3, xlab = "Distance", ylab = "Variance",
+ ylim = c(0, ymax), cex.lab = 1.2)
+ lines(vario$Dist, vario$All + z * vario$se, lty = 1)
+ lines(vario$Dist, vario$All - z * vario$se, lty = 1)
+ lines(vario$Dist, vario$Sum, type = "b", lty = 2,
+ pch = 3)
+ for (i in 6:(b + 3)) {
+ lines(vario$Dist, vario[, i], type = "b", lty = 1,
+ pch = i - 6)
+ points(x = 1.2 * grain,
+ y = ymax - ymax * (b + 6 - i)/20, pch = i - 6)
+ }
+ text(x = rep(2 * grain, b - 1), y = ymax - ymax *
+ c(2:b)/20, label = label[c(2, b:3) + 3], pos = 4,
+ cex = 1.2)
+ points(x = 1.2 * grain, y = ymax - ymax * 2/20, pch = 3)
+ for (i in 2:b) {
+ lines(x = c(0.7, 1.1) * grain,
+ y = rep(ymax - ymax * i/20, 2),
+ lty = c(1, 2, 1, 1, 1)[i])
+ lines(x = c(1.3, 1.7) * grain,
+ y = rep(ymax - ymax * i/20, 2), lty = c(1, 2, 1, 1, 1)[i])
+ }
+ text(x = c(vario$Dist), y = rep(0, length(vario$Dist)),
+ label = c(vario$n), cex = 0.8)
+ lines(x = rep(max(object$H)/2, 2), y = c(-10, ymax +
+ 10), lty = 3)
+ text(x = 2 * grain, y = ymax - ymax/20, label = "C.I. for total variance",
+ pos = 4, cex = 1.2)
+ lines(x = c(0.7, 1.7) * grain, y = rep(ymax - ymax/20, 2),
+ lty = 1)
+ }
+ else {
+ plot(vario$Dist, vario$All, type = "b", lty = 1,
+ pch = 0, xlab = "Distance", ylab = "Variance",
+ ylim = c(0, ymax), cex.lab = 1.2)
+ lines(c(0, 10), rep(object$tot.chi, 2), lty = 5)
+ lines(x = c(0.7, 1.7) * grain, y = rep(ymax - ymax *
+ b/20, 2), lty = 5)
+ text(x = 2 * grain, y = ymax - ymax * b/20, label = "Global variance estimate",
+ pos = 4, cex = 1.2)
+ text(x = c(vario$Dist), y = rep(0, length(vario$Dist)),
+ label = c(vario$n), cex = 0.8)
+ lines(x = rep(max(object$H)/2, 2), y = c(-10, ymax +
+ 10), lty = 3)
+ text(x = 2 * grain, y = ymax - ymax/20, label = "Total variance",
+ pos = 4, cex = 1.2)
+ lines(x = c(0.7, 1.7) * grain, y = rep(ymax - ymax/20,
+ 2), lty = 1)
+ }
+ }
+ if (is.numeric(object$vario$CA.signif)) {
+ a <- c(1:nrow(object$vario))[object$vario$CA.signif <
+ alpha]
+ points(vario$Dist[a], object$vario$CA[a], pch = 15)
+ points(x = 1.2 * grain, y = ymax - ymax * (b + 1)/20,
+ pch = 15)
+ text(x = 2 * grain, y = ymax - ymax * (b + 1)/20, pos = 4,
+ cex = 1.2, label = c("Sign. autocorrelation"))
+ if (is.numeric(object$CCA$rank)) {
+ inflation <- 1 - weighted.mean(object$vario$CA, object$vario$n)/
+ weighted.mean(object$vario$CA[-a],
+ object$vario$n[-a])
+ cat("Error variance of regression model underestimated by",
+ round(inflation * 100, 1), "percent", "\n")
+ }
+ }
+ invisible()
+}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-05-07 13:30:37 UTC (rev 338)
+++ pkg/inst/ChangeLog 2008-05-07 13:54:12 UTC (rev 339)
@@ -4,6 +4,11 @@
Version 1.12-14 (opened May 7, 2008)
+ * mso: added Helene Wagner's mso function for direct multi-scale
+ ordination with constrained ordination or spatial partioning of
+ 'cca' and 'rda' results (Ecology 85, 342-351; 2004). Thanks to
+ Helene Wagner for allowing the inclusion of the code.
+
* MOStest: new function to implement Mitchell-Olds & Shaw test for
the location of quadratic extreme in a defined interval.
Added: pkg/man/mso.Rd
===================================================================
--- pkg/man/mso.Rd (rev 0)
+++ pkg/man/mso.Rd 2008-05-07 13:54:12 UTC (rev 339)
@@ -0,0 +1,115 @@
+\name{mso}
+\alias{mso}
+\alias{plot.mso}
+
+\title{ Functions for performing and displaying a spatial partitioning
+of cca or rda results}
+
+\description{ The function \code{mso} adds an attribute \code{vario} to
+ an object of class \code{"cca"} that describes the spatial
+ partitioning of the \code{\link{cca}} object and performs an optional
+ permutation test for the spatial independence of residuals. The
+ function \code{plot.mso} creates a diagnostic plot of the spatial
+ partitioning of the \code{"cca"} object. }
+
+\usage{
+mso(object.cca, object.xy, grain = 1, round.up = FALSE, permutations = FALSE)
+\method{plot}{mso}(x, alpha = 0.05, explained = FALSE, ...)
+}
+\arguments{
+ \item{object.cca}{ An object of class cca, created by the \code{\link{cca}} or
+ \code{\link{rda}} function.}
+ \item{object.xy}{ A vector, matrix or data frame with the spatial
+ coordinates of the data represented by object.cca. Must have the
+ same number of rows as \code{object.cca\$CA\$Xbar} (see
+ \code{\link{cca.object}}). }
+ \item{grain}{ Interval size for distance classes.}
+ \item{round.up}{ Determines the choice of breaks. If false, distances
+ are rounded to the nearest multiple of grain. If true, distances are
+ rounded to the upper multiple of grain.}
+ \item{permutations}{ If false, suppresses the permutation test. If an
+ integer, determines the number of permutations for the Mantel test
+ of spatial independence of residual inertia.}
+ \item{x}{A result object of \code{mso}.}
+ \item{alpha}{ Significance level for the two-sided permutation test of
+ the Mantel statistic for spatial independence of residual inertia
+ and for the point-wise envelope of the variogram of the total
+ variance. A Bonferroni-type correction can be achieved by dividing
+ the overall significance value (e.g. 0.05) by the number of distance
+ classes.}
+ \item{explained}{ If false, suppresses the plotting of the variogram
+ of explained variance.}
+ \item{\dots}{Other arguments passed to functions.}
+}
+\details{
+ The Mantel test is an adaptation of the function \code{\link{mantel}} of the
+ \pkg{vegan} package to the parallel testing of several distance classes. It
+ compares the mean inertia in each distance class to the pooled mean
+ inertia of all other distance classes.
+
+ If there are explanatory variables (RDA, CCA, pRDA, pCCA) and a
+ significance test for residual autocorrelation was performed when
+ running the function \code{mso}, the function \code{plot.mso} will
+ print an estimate of how much the autocorrelation (based on
+ significant distance classes) causes the global error variance of the
+ regression analysis to be underestimated
+
+}
+\value{
+ The function \code{mso} returns an amended \code{cca} or \code{rda}
+ object with the additional attributes \code{grain}, \code{H},
+ \code{H.test} and \code{vario}.
+ \item{grain }{ The grain attribute defines the interval size of the
+ distance classes .}
+ \item{H }{ H is an object of class 'dist' and contains the geographic
+ distances between observations.}
+ \item{H.test }{ H.test contains a set of dummy variables that describe
+ which pairs of observations (rows = elements of \code{object\$H}) fall in
+ which distance class (columns). }
+ \item{vario }{ The vario attribute is a data frame that contains some
+ or all of the following components for the rda case (cca case in
+ brackets):
+ \item{H }{ Distance class as multiples of grain.}
+ \item{Dist }{ Average distance of pairs of observations in distance class H.}
+ \item{n }{ Number of unique pairs of observations in distance class
+ H.}
+ \item{All }{ Empirical (chi-square) variogram of total variance
+ (inertia).}
+ \item{Sum }{ Sum of empirical (chi-square) variograms of explained
+ and residual variance (inertia).}
+ \item{CA }{ Empirical (chi-square) variogram of residual variance
+ (inertia).}
+ \item{CCA }{ Empirical (chi-square) variogram of explained variance
+ (inertia).}
+ \item{pCCA }{ Empirical (chi-square) variogram of conditioned
+ variance (inertia).}
+ \item{se }{ Standard error of the empirical (chi-square) variogram
+ of total variance (inertia).}
+ \item{CA.signif }{P-value of permutation test for spatial
+ independence of residual variance (inertia).}
+ }
+}
+\references{ Wagner, H.H. 2004. Direct multi-scale ordination with
+ canonical correspondence analysis. \emph{Ecology} 85: 342--351. }
+\author{ The responsible author was Helene Wagner.}
+\note{ The function is based on the code published in the Ecological
+ Archives E085-006 (\url{http://www.esapubs.org/archive/ecol/E085/006/default.htm}). }
+\seealso{ Function \code{\link{cca}} and \code{\link{rda}},
+ \code{\link{cca.object}}. }
+\examples{
+## Reconstruct worked example of Wagner (submitted):
+X <- matrix(c(1, 2, 3, 2, 1, 0), 3, 2)
+Y <- c(3, -1, -2)
+tmat <- c(1:3)
+## Canonical correspondence analysis (cca):
+Example.cca <- cca(X, Y)
+Example.cca <- mso(Example.cca, tmat)
+plot(Example.cca)
+Example.cca$vario
+
+## Correspondence analysis (ca):
+Example.ca <- mso(cca(X), tmat)
+plot(Example.ca)
+}
+\keyword{ spatial }
+\keyword{ multivariate }
More information about the Vegan-commits
mailing list