[Vegan-commits] r391 - in branches/1.13: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 4 08:54:51 CEST 2008
Author: jarioksa
Date: 2008-06-04 08:54:50 +0200 (Wed, 04 Jun 2008)
New Revision: 391
Added:
branches/1.13/R/msoplot.R
branches/1.13/R/print.mso.R
Removed:
branches/1.13/R/plot.mso.R
Modified:
branches/1.13/R/mso.R
branches/1.13/inst/ChangeLog
branches/1.13/man/mso.Rd
Log:
mso in branches/1.13 updated to r389 in pkg/
Modified: branches/1.13/R/mso.R
===================================================================
--- branches/1.13/R/mso.R 2008-06-04 06:45:21 UTC (rev 390)
+++ branches/1.13/R/mso.R 2008-06-04 06:54:50 UTC (rev 391)
@@ -1,17 +1,23 @@
-mso <- function (object.cca, object.xy, grain = 1, round.up = FALSE,
- permutations = FALSE)
+`mso` <-
+ function (object.cca, object.xy, grain = 1, round.up = FALSE,
+ permutations = FALSE)
{
+ if (inherits(object.cca, "mso")) {
+ rm <- which(class(object.cca) == "mso")
+ class(object.cca) <- class(object.cca)[-rm]
+ }
object <- object.cca
xy <- object.xy
N <- nrow(object$CA$Xbar)
- if (class(object) == "rda")
+ if (inherits(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
+ hmax <- round((max(Dist)/2)/grain) *grain
+ H[H > hmax] <- max(H)
object$H <- H
H <- as.vector(H)
Dist <- sapply(split(Dist, H), mean)
@@ -71,6 +77,8 @@
object$vario$CA.signif <- apply((perm >= matrix(statistic,
nrow(perm), ncol(perm)))/permutations, 1, sum)
}
+ object$call <- match.call()
class(object) <- c("mso", class(object))
object
}
+
Copied: branches/1.13/R/msoplot.R (from rev 390, pkg/R/msoplot.R)
===================================================================
--- branches/1.13/R/msoplot.R (rev 0)
+++ branches/1.13/R/msoplot.R 2008-06-04 06:54:50 UTC (rev 391)
@@ -0,0 +1,97 @@
+`msoplot` <-
+ 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")
+ ## You should not change par, or at least you must put
+ ## back the old values when exiting:
+ ## op <- par(omi = c(0.5, 0.5, 0, 0))
+ ## on.exit(par(op))
+ ##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()
+}
+
Deleted: branches/1.13/R/plot.mso.R
===================================================================
--- branches/1.13/R/plot.mso.R 2008-06-04 06:45:21 UTC (rev 390)
+++ branches/1.13/R/plot.mso.R 2008-06-04 06:54:50 UTC (rev 391)
@@ -1,91 +0,0 @@
-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()
-}
Copied: branches/1.13/R/print.mso.R (from rev 390, pkg/R/print.mso.R)
===================================================================
--- branches/1.13/R/print.mso.R (rev 0)
+++ branches/1.13/R/print.mso.R 2008-06-04 06:54:50 UTC (rev 391)
@@ -0,0 +1,9 @@
+`print.mso` <-
+ function(x, digits = max(3, getOption("digits") - 3), ...)
+{
+ NextMethod(x, "print", digits = digits, ...)
+ cat("mso variogram:\n\n")
+ print(x$vario, digits = digits, ...)
+ invisible(x)
+}
+
Modified: branches/1.13/inst/ChangeLog
===================================================================
--- branches/1.13/inst/ChangeLog 2008-06-04 06:45:21 UTC (rev 390)
+++ branches/1.13/inst/ChangeLog 2008-06-04 06:54:50 UTC (rev 391)
@@ -5,6 +5,9 @@
Version 1.13-1 (opened May 22, 2008)
+ * merged mso revisions since r372: now similar as in devel branch
+ at r389.
+
* merged 371: mite.xy data.
* merged 368: method name to taxa2dist
Modified: branches/1.13/man/mso.Rd
===================================================================
--- branches/1.13/man/mso.Rd 2008-06-04 06:45:21 UTC (rev 390)
+++ branches/1.13/man/mso.Rd 2008-06-04 06:54:50 UTC (rev 391)
@@ -1,6 +1,7 @@
\name{mso}
\alias{mso}
-\alias{plot.mso}
+\alias{print.mso}
+\alias{msoplot}
\title{ Functions for performing and displaying a spatial partitioning
of cca or rda results}
@@ -14,7 +15,7 @@
\usage{
mso(object.cca, object.xy, grain = 1, round.up = FALSE, permutations = FALSE)
-\method{plot}{mso}(x, alpha = 0.05, explained = FALSE, ...)
+msoplot(x, alpha = 0.05, explained = FALSE, ...)
}
\arguments{
\item{object.cca}{ An object of class cca, created by the \code{\link{cca}} or
@@ -104,12 +105,30 @@
## Canonical correspondence analysis (cca):
Example.cca <- cca(X, Y)
Example.cca <- mso(Example.cca, tmat)
-plot(Example.cca)
+msoplot(Example.cca)
Example.cca$vario
## Correspondence analysis (ca):
Example.ca <- mso(cca(X), tmat)
-plot(Example.ca)
+msoplot(Example.ca)
+
+## Unconstrained ordination with test for autocorrelation
+## using oribatid mite data set as in Wagner (2004)
+data(mite)
+data(mite.env)
+data(mite.xy)
+
+mite.cca <- cca(log(mite + 1))
+mite.cca <- mso(mite.cca, mite.xy, grain = 1, permutations = 100)
+msoplot(mite.cca)
+mite.cca
+
+## Constrained ordination with test for residual autocorrelation
+## and scale-invariance of species-environment relationships
+mite.cca <- cca(log(mite + 1) ~ SubsDens + WatrCont + Substrate + Shrub + Topo, mite.env)
+mite.cca <- mso(mite.cca, mite.xy, permutations = 100)
+msoplot(mite.cca)
+mite.cca
}
\keyword{ spatial }
\keyword{ multivariate }
More information about the Vegan-commits
mailing list