[Vegan-commits] r1086 - in pkg/vegan: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Dec 19 16:14:38 CET 2009
Author: gsimpson
Date: 2009-12-19 16:14:37 +0100 (Sat, 19 Dec 2009)
New Revision: 1086
Added:
pkg/vegan/R/ordimedian.R
Modified:
pkg/vegan/R/betadisper.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/betadisper.Rd
Log:
Added spatial median option to betadisper, plus associated code.
Modified: pkg/vegan/R/betadisper.R
===================================================================
--- pkg/vegan/R/betadisper.R 2009-12-16 16:59:29 UTC (rev 1085)
+++ pkg/vegan/R/betadisper.R 2009-12-19 15:14:37 UTC (rev 1086)
@@ -7,8 +7,6 @@
if(!inherits(d, "dist"))
stop("distances 'd' must be a 'dist' object")
type <- match.arg(type)
- if(type == "median")
- .NotYetUsed('type = "median"')
## checks for groups - need to be a factor for later
if(!is.factor(group))
group <- as.factor(group)
@@ -50,7 +48,11 @@
## store which are the positive eigenvalues
pos <- eig > 0
## group centroids in PCoA space
- centroids <- apply(vectors, 2, function(x) tapply(x, group, mean))
+ centroids <-
+ switch(type,
+ centroid = apply(vectors, 2, function(x) tapply(x, group, mean)),
+ median = ordimedian(vectors, group, choices = 1:ncol(vectors))
+ )
## for each of the groups, calculate distance to centroid for
## observation in the group
if(is.matrix(centroids)) {
@@ -65,8 +67,7 @@
dist.neg <- 0
}
} else {
- dist.pos <- sweep(vectors[, pos, drop=FALSE], 2,
- centroids[pos])
+ dist.pos <- sweep(vectors[, pos, drop=FALSE], 2, centroids[pos])
dist.pos <- rowSums(dist.pos^2)
if (any(!pos)) {
dist.neg <- sweep(vectors[, !pos, drop=FALSE], 2,
@@ -89,5 +90,6 @@
group = group, centroids = centroids, call = match.call())
class(retval) <- "betadisper"
attr(retval, "method") <- attr(d, "method")
+ attr(retval, "type") <- type
retval
}
Added: pkg/vegan/R/ordimedian.R
===================================================================
--- pkg/vegan/R/ordimedian.R (rev 0)
+++ pkg/vegan/R/ordimedian.R 2009-12-19 15:14:37 UTC (rev 1086)
@@ -0,0 +1,37 @@
+## Ordimedian finds the spatial medians for groups. Spatial medians
+## are L1 norms or statistics that minimize sum of distances of points
+## from the statistic and 1d they are the medians. The current
+## algorithm minimizes the L1 norm with optim and is pretty
+## inefficient. Package ICSNP has a better algorithm (and we may steal
+## it from them later).
+`ordimedian` <-
+ function(ord, groups, display = "sites", label = FALSE, ...)
+{
+ ## Sum of distances from the statistic
+ medfun <-
+ function(x, ord) sum(sqrt(rowSums(sweep(ord, 2, x)^2)),
+ na.rm = TRUE)
+ ## derivative of medfun (if NULL, optim will use numerical
+ ## differentiation)
+ dmedfun <- function(x, ord) {
+ up <- -sweep(ord, 2, x)
+ dn <- sqrt(rowSums(sweep(ord, 2, x)^2))
+ colSums(sweep(up, 1, dn, "/"))
+ }
+ #dmedfun <- NULL
+ pts <- scores(ord, display = display, ...)
+ inds <- names(table(groups))
+ medians <- matrix(0, nrow = length(inds), ncol = ncol(pts))
+ rownames(medians) <- inds
+ colnames(medians) <- colnames(pts)
+ for (i in inds) {
+ X <- pts[groups == i, , drop = FALSE]
+ medians[i, ] <- optim(apply(X, 2, median, na.rm = TRUE),
+ fn = medfun, gr = dmedfun,
+ ord = X, method = "BFGS")$par
+ if(label)
+ ordiArgAbsorber(medians[i,1], medians[i,2], label = i,
+ FUN = text, ...)
+ }
+ invisible(medians)
+}
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2009-12-16 16:59:29 UTC (rev 1085)
+++ pkg/vegan/inst/ChangeLog 2009-12-19 15:14:37 UTC (rev 1086)
@@ -4,14 +4,14 @@
Version 1.16-34 (opened December 13, 2009) -- (codename: Lucia)
- * twostagechao: deleted, because question regarding the
+ * twostagechao: deleted, because question regarding the
notation in the original article remained unanswered.
* permat* cleanup: permat.control deleted, not needed any longer due to changes
in oecosimu, documentation for plot.permat merged with permat* help page.
* goodness.cca: was not using correct weights. Broken in r992 (Sep
- 6, 2009).
+ 6, 2009).
* vegdist: added UI for the alternative Gower index ("altGower")
that skips double zeros. The C code has been in vegan since Oct 24
@@ -26,14 +26,19 @@
places in permutest.cca and in anova.ccabyterm.
* permutest.cca: uses pasteCall() in print.
-
+
+ * betadisper: added spatial median code written by Jari to betadisper.
+ All associated methods work OK, and this should be the default method
+ but needs a little checking first. This requires the new function
+ 'ordimedian', which is currently undocumented in ?betadisper.
+
Version 1.16-33 (closed December 13, 2009)
- * DESCRIPTION: Guillaume Blanchet added to the Authors.
+ * DESCRIPTION: Guillaume Blanchet added to the Authors.
* spandepth: New function to find the depths of nodes in a
spanning tree from spantree() function. Per request of
- W. E. Sharp.
+ W. E. Sharp.
* bstick: refuses to analyse capscale() models with imaginary
components because I have no idea how to do it. Reported by
Modified: pkg/vegan/man/betadisper.Rd
===================================================================
--- pkg/vegan/man/betadisper.Rd 2009-12-16 16:59:29 UTC (rev 1085)
+++ pkg/vegan/man/betadisper.Rd 2009-12-19 15:14:37 UTC (rev 1086)
@@ -6,6 +6,7 @@
\alias{boxplot.betadisper}
\alias{print.betadisper}
\alias{TukeyHSD.betadisper}
+\alias{ordimedian}
\title{Multivariate homogeneity of groups dispersions (variances)}
\description{
@@ -25,7 +26,7 @@
Tukey's 'Honest Significant Difference' method.
}
\usage{
-betadisper(d, group, type = c("centroid", "median"))
+betadisper(d, group, type = c("centroid","median"))
\method{anova}{betadisper}(object, \dots)
@@ -49,8 +50,8 @@
or an object that can be coerced to a factor using
\code{\link[base]{as.factor}}. Can consist of a factor with a single
level (i.e.~one group).}
- \item{type}{the type of analysis to perform. Only \code{type =
- "centroid"} is currently supported.}
+ \item{type}{the type of analysis to perform. Use the spatial median or
+ the group centroid?}
\item{display}{character; partial match to access scores for
\code{"sites"} or \code{"species"}.}
\item{object, x}{an object of class \code{"betadisper"}, the result of a
@@ -222,11 +223,20 @@
dis[c(2, 20)] <- NA
mod2 <- betadisper(dis, groups) ## warnings
mod2
-permutest(mod, control = permControl(nperm = 100))
+permutest(mod2, control = permControl(nperm = 100))
anova(mod2)
plot(mod2)
boxplot(mod2)
plot(TukeyHSD(mod2))
+
+## Using spatial median
+mod3 <- betadisper(dis, groups, type = "median")
+mod3
+permutest(mod3, control = permControl(nperm = 100))
+anova(mod3)
+plot(mod23)
+boxplot(mod3)
+plot(TukeyHSD(mod3))
}
\keyword{methods}
\keyword{multivariate}
More information about the Vegan-commits
mailing list