[Vegan-commits] r257 - in pkg: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 9 14:29:40 CET 2008
Author: gsimpson
Date: 2008-03-09 14:29:39 +0100 (Sun, 09 Mar 2008)
New Revision: 257
Added:
pkg/R/TukeyHSD.betadisper.R
pkg/man/TukeyHSD.betadisper.Rd
pkg/man/permDisper.Rd
Modified:
pkg/DESCRIPTION
pkg/R/betadisper.R
pkg/R/permDisper.R
pkg/R/print.permDisper.R
pkg/inst/ChangeLog
pkg/man/betadisper.Rd
Log:
Adds pairwise comparisons in permDisper() and TukeyHSD.betadisper()
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/DESCRIPTION 2008-03-09 13:29:39 UTC (rev 257)
@@ -1,7 +1,7 @@
Package: vegan
Title: Community Ecology Package
-Version: 1.12-3
-Date: Mar 7, 2008
+Version: 1.12-4
+Date: Mar 9, 2008
Author: Jari Oksanen, Roeland Kindt, Pierre Legendre, Bob O'Hara, Gavin L. Simpson,
M. Henry H. Stevens
Maintainer: Jari Oksanen <jari.oksanen at oulu.fi>
Added: pkg/R/TukeyHSD.betadisper.R
===================================================================
--- pkg/R/TukeyHSD.betadisper.R (rev 0)
+++ pkg/R/TukeyHSD.betadisper.R 2008-03-09 13:29:39 UTC (rev 257)
@@ -0,0 +1,7 @@
+`TukeyHSD.betadisper` <- function(x, which = "group", ordered = FALSE,
+ conf.level = 0.95, ...) {
+ df <- data.frame(distances = x$distances, group = x$group)
+ mod.aov <- aov(distances ~ group, data = df)
+ TukeyHSD(mod.aov, which = which, ordered = ordered,
+ conf.level = conf.level, ...)
+}
Modified: pkg/R/betadisper.R
===================================================================
--- pkg/R/betadisper.R 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/R/betadisper.R 2008-03-09 13:29:39 UTC (rev 257)
@@ -23,11 +23,11 @@
centroids <- apply(vectors, 2, function(x) tapply(x, group, mean))
## for each of the groups, calculate distance to centroid for
## observation in the group
- dist.pos <- vectors[, pos, drop=FALSE] - centroids[group, pos, drop = FALSE]
+ dist.pos <- vectors[, pos, drop=FALSE] - centroids[group, pos, drop=FALSE]
dist.pos <- rowSums(dist.pos^2)
if (any(!pos)) {
dist.neg <- vectors[, !pos, drop=FALSE] -
- centroids[group, !pos, drop = FALSE]
+ centroids[group, !pos, drop=FALSE]
dist.neg <- rowSums(dist.neg^2)
} else {
dist.neg <- 0
@@ -35,7 +35,7 @@
## zij are the distances of each point to its group centroid
zij <- sqrt(dist.pos - dist.neg)
## add in correct labels
- colnames(vectors) <- colnames(centroids) <- names(eig) <-
+ colnames(vectors) <- colnames(centroids) <- names(eig) <-
paste("PCoA", 1:n, sep = "")
rownames(vectors) <- names(zij) <- attr(d, "Labels")
retval <- list(eig = eig, vectors = vectors, distances = zij,
Modified: pkg/R/permDisper.R
===================================================================
--- pkg/R/permDisper.R 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/R/permDisper.R 2008-03-09 13:29:39 UTC (rev 257)
@@ -1,5 +1,18 @@
-`permDisper` <- function(object, control = permControl(nperm = 999))
+`permDisper` <- function(object, pairwise = FALSE,
+ control = permControl(nperm = 999))
{
+ t.statistic <- function(x, y) {
+ m <- length(x)
+ n <- length(y)
+ xbar <- sum(x) / m
+ ybar <- sum(y) / n
+ #xvar <- var(x)
+ #yvar <- var(y)
+ xvar <- .Internal(cov(x, NULL, 1, FALSE))
+ yvar <- .Internal(cov(y, NULL, 1, FALSE))
+ pooled <- sqrt(((m-1)*xvar + (n-1)*yvar) / (m+n-2))
+ (xbar - ybar) / (pooled * sqrt(1/m + 1/n))
+ }
if(!inherits(object, "betadisper"))
stop("Only for class \"betadisper\"")
nobs <- length(object$distances)
@@ -9,8 +22,19 @@
resids <- qr.resid(mod.Q, object$distances)
res <- numeric(length = control$nperm + 1)
res[1] <- summary(mod)$fstatistic[1]
+ ## pairwise comparisons
+ if(pairwise) {
+ ## unique pairings
+ combin <- combn(levels(object$group), 2)
+ n.pairs <- ncol(combin)
+ t.stats <- matrix(0, ncol = n.pairs, nrow = control$nperm + 1)
+ t.stats[1,] <- apply(combn(levels(object$group), 2), 2, function(x) {
+ t.statistic(object$distances[object$group == x[1]],
+ object$distances[object$group == x[2]])})
+ }
for(i in seq(along = res[-1])) {
- perm.resid <- resids[permuted.index2(nobs, control = control)]
+ perm <- permuted.index2(nobs, control = control)
+ perm.resid <- resids[perm]
f <- qr.fitted(mod.Q, perm.resid)
mss <- sum((f - mean(f))^2)
r <- qr.resid(mod.Q, perm.resid)
@@ -18,14 +42,35 @@
rdf <- nobs - p
resvar <- rss / rdf
res[i+1] <- (mss / (p - 1)) / resvar
+ ## pairwise comparisons
+ if(pairwise) {
+ for(j in seq_len(n.pairs)) {
+ grp1 <- object$distance[perm][object$group == combin[1, j]]
+ grp2 <- object$distance[perm][object$group == combin[2, j]]
+ t.stats[i+1, j] <- t.statistic(grp1, grp2)
+ }
+ }
}
pval <- sum(res >= res[1]) / length(res)
+ if(pairwise) {
+ df <- apply(combin, 2, function(x) {
+ length(object$distances[object$group == x[1]]) +
+ length(object$distance[object$group == x[2]]) - 2})
+ pairwise <- list(observed = 2 * pt(-abs(t.stats[1,]), df),
+ permuted = apply(t.stats, 2,
+ function(x) sum(abs(x) >= x[1])/length(x)))
+ names(pairwise$observed) <- names(pairwise$permuted) <-
+ apply(combin, 2, paste, collapse = "-")
+ } else {
+ pairwise <- NULL
+ }
mod.aov <- anova(object)
retval <- cbind(mod.aov[, 1:4], c(control$nperm, NA), c(pval, NA))
dimnames(retval) <- list(c("Groups", "Residuals"),
c("Df", "Sum Sq", "Mean Sq", "F", "N.Perm",
"Pr(>F)"))
- retval <- list(tab = retval, control = control)
- class(retval) <- "permDisper" ##c("permDisper", class(retval))
+ retval <- list(tab = retval, pairwise = pairwise,
+ groups = levels(object$group), control = control)
+ class(retval) <- "permDisper"
retval
}
Modified: pkg/R/print.permDisper.R
===================================================================
--- pkg/R/print.permDisper.R 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/R/print.permDisper.R 2008-03-09 13:29:39 UTC (rev 257)
@@ -26,5 +26,15 @@
signif.stars = getOption("show.signif.stars"),
has.Pvalue = has.P, P.values = has.P, cs.ind = NULL,
zap.ind = zap.i, tst.ind = tst.i, na.print = "", ...)
+ if(!is.null(x$pairwise)) {
+ cat("\nPairwise comparisons:", sep = "\n")
+ writeLines(strwrap("(Observed p-value below diagonal,\npermuted p-value above diagonal)\n"))
+ n.grp <- length(x$groups)
+ mat <- matrix(NA, ncol = n.grp, nrow = n.grp)
+ colnames(mat) <- rownames(mat) <- x$groups
+ mat[lower.tri(mat)] <- x$pairwise$observed
+ mat[upper.tri(mat)] <- x$pairwise$permuted
+ printCoefmat(mat, na.print = "", digits = digits)
+ }
invisible(x)
}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/inst/ChangeLog 2008-03-09 13:29:39 UTC (rev 257)
@@ -2,11 +2,22 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 1.2-3 (Mar 8, 2008, working...)
+Version 1.12-4 (opened Mar 9, 2008)
+
+ * permDisper: Added pairwise comparisons of group dispersions
+ via a classical t test and via permutation test, accessed via
+ new argument 'pairwise = TRUE'. 'permDisper' is now documented
+ seperately from 'betadisper'.
+ * TukeyHSD.betadisper: A method for 'TukeyHSD' to calculate
+ Tukey's Honest Significant Differences for the grouping factor
+ in 'betadisper'.
+
+Version 1.12-3 (closed Mar 9, 2008)
+
* betadiver: a new function that implements all indices of beta
diversity reviewed by Koleff et al. (J. Anim. Ecol., 72, 367-382;
- 2003), with a plot function to produce triangular plots.
+ 2003), with a plot function to produce triangular plots.
* isomap: added dynamic, spinnable 3D graphics using rgl
(rgl.isomap).
Added: pkg/man/TukeyHSD.betadisper.Rd
===================================================================
--- pkg/man/TukeyHSD.betadisper.Rd (rev 0)
+++ pkg/man/TukeyHSD.betadisper.Rd 2008-03-09 13:29:39 UTC (rev 257)
@@ -0,0 +1,50 @@
+\name{TukeyHSD.betadisper}
+\alias{TukeyHSD.betadisper}
+\title{Compute Tukey Honest Significant Differences}
+
+\description{
+ Create a set of confidence intervals on the differences between
+ the means of the levels of the grouping factor with the specified
+ family-wise probability of coverage. The intervals are based on the
+ Studentized range statistic, Tukey's 'Honest Significant
+ Difference' method.
+}
+\usage{
+\method{TukeyHSD}{betadisper}(x, which = "group", ordered = FALSE,
+ conf.level = 0.95, \ldots)
+}
+\arguments{
+ \item{x}{A fitted model object of class \code{\link{betadisper}}.}
+ \item{which}{A character vector listing terms in the fitted model for
+ which the intervals should be calculated. Defaults to the grouping
+ factor.}
+ \item{ordered}{Logical; see \code{\link{TukeyHSD}}.}
+ \item{conf.level}{A numeric value between zero and one giving the
+ family-wise confidence level to use.}
+ \item{\ldots}{Additional arguments passed to \code{\link{TukeyHSD}}.}
+}
+\details{
+ This is a simple wrapper to \code{\link{TukeyHSD.aov}}. The user is
+ directed to read the help file for \code{\link{TukeyHSD}} before using
+ this function.
+
+ In particular, note the statement about using the function with
+ unbalanced designs.
+}
+\author{
+ Douglas Bates
+}
+\seealso{
+ \code{\link{betadisper}}, \code{\link{permDisper}},
+ \code{\link{aov}}.
+}
+\examples{
+## continue the example from ?betadisper
+example(betadisper)
+
+mod.hsd <- TukeyHSD(mod)
+plot(mod.hsd)
+
+}
+\keyword{models}
+\keyword{design}
\ No newline at end of file
Modified: pkg/man/betadisper.Rd
===================================================================
--- pkg/man/betadisper.Rd 2008-03-08 16:20:20 UTC (rev 256)
+++ pkg/man/betadisper.Rd 2008-03-09 13:29:39 UTC (rev 257)
@@ -1,12 +1,10 @@
\name{betadisper}
\alias{betadisper}
-\alias{permDisper}
\alias{scores.betadisper}
\alias{anova.betadisper}
\alias{plot.betadisper}
\alias{boxplot.betadisper}
\alias{print.betadisper}
-\alias{print.permDisper}
\title{Multivariate homogeneity of groups dispersions (variances)}
\description{
Implements Marti Anderson's PERMDISP2 procedure for the analysis of
@@ -22,8 +20,6 @@
\method{anova}{betadisper}(object, \dots)
-permDisper(object, control = permControl(nperm = 999))
-
\method{scores}{betadisper}(x, display = c("sites", "centroids"),
choices = c(1,2), \dots)
@@ -42,13 +38,10 @@
\code{\link[base]{as.factor}}.}
\item{type}{the type of analysis to perform. Only \code{type =
"centroid"} is currently supported.}
+ \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
call to \code{betadisper}.}
- \item{control}{a list of control values for the permutations
- to replace the default values returned by the function
- \code{\link{permControl}}.}
- \item{display}{character; partial match to access scores for
- \code{"sites"} or \code{"species"}.}
\item{choices, axes}{the principal coordinate axes wanted.}
\item{hull}{logical; should the convex hull for each group be plotted?}
\item{cex, ylab, xlab, main, sub}{graphical parameters. For details,
@@ -83,10 +76,16 @@
To test if one or more groups is more variable than the others, ANOVA
of the distances to group centroids can be performed and parametric
theory used to interpret the significance of F. An alternative is to
- use a permutation test. \code{permDisper} permutes model residuals to
- generate a permutation distribution of F under the Null hypothesis of
- no difference in dispersion between groups.
+ use a permutation test. \code{\link{permDisper}} permutes model
+ residuals to generate a permutation distribution of F under the Null
+ hypothesis of no difference in dispersion between groups.
+ Pairwise comprisons of group mean dispersions can also be performed
+ using \code{\link{permDisper}}. An alternative to the classical
+ comparison of group dispersions, is to calculate Tukey's Honest
+ Significant Differences between groups, via
+ \code{\link{TukeyHSD.betadisper}}.
+
The results of the analysis can be visualised using the \code{plot}
and \code{boxplot} methods.
@@ -98,11 +97,6 @@
The \code{anova} method returns an object of class \code{"anova"}
inheriting from class \code{"data.frame"}.
- \code{permDisper} returns an object of class \code{"permDisper"}, a list
- with components \code{tab}, the ANOVA table which is an object
- inheriting from class \code{"data.frame"}, and \code{control}, the
- result of a call to \code{\link{permControl}}..
-
The \code{scores} method returns a list with one or both of the
components \code{"sites"} and \code{"centroids"}.
@@ -138,8 +132,10 @@
\strong{9(6)}, 683--693.
}
\author{Gavin L. Simpson}
-\seealso{\code{\link[stats]{anova.lm}}, \code{\link{scores}},
- \code{\link[graphics]{boxplot}}, \code{\link{betadiver}}.}
+\seealso{\code{\link{permDisper}}, \code{\link[stats]{anova.lm}},
+ \code{\link{scores}}, \code{\link[graphics]{boxplot}},
+ \code{\link{TukeyHSD.betadisper}}. Further measure of beta diversity
+ can be found in \code{\link{betadiver}}.}
\examples{
data(varespec)
@@ -157,8 +153,12 @@
anova(mod)
## Permutation test for F
-permDisper(mod)
+permDisper(mod, pairwise = TRUE)
+## Tukey's Honest Significant Differences
+(mod.HSD <- TukeyHSD(mod))
+plot(mod.HSD)
+
## Plot the groups and distances to centroids on the
## first two PCoA axes
plot(mod)
Added: pkg/man/permDisper.Rd
===================================================================
--- pkg/man/permDisper.Rd (rev 0)
+++ pkg/man/permDisper.Rd 2008-03-09 13:29:39 UTC (rev 257)
@@ -0,0 +1,88 @@
+\name{permDisper}
+\alias{permDisper}
+\alias{print.permDisper}
+\title{Permutation test of ultivariate homogeneity of groups dispersions
+ (variances)}
+\description{
+ Implements a permutation-based test of multivariate homogeneity of
+ group dispersions (variances) for the results of a call to
+ \code{\link{betadisper}}.
+}
+\usage{
+permDisper(object, pairwise = FALSE, control = permControl(nperm = 999))
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{object}{an object of class \code{"betadisper"}, the result of a
+ call to \code{betadisper}.}
+ \item{pairwise}{logical; perform pairwise comparisons of group means?}
+ \item{control}{a list of control values for the permutations
+ to replace the default values returned by the function
+ \code{\link{permControl}}.}
+}
+\details{
+ To test if one or more groups is more variable than the others, ANOVA
+ of the distances to group centroids can be performed and parametric
+ theory used to interpret the significance of F. An alternative is to
+ use a permutation test. \code{permDisper} permutes model residuals to
+ generate a permutation distribution of F under the Null hypothesis of
+ no difference in dispersion between groups.
+
+ Pairwise comprisons of group mean dispersions can be performed by
+ setting argument \code{pairwise} to \code{TRUE}. A classicial t test
+ is performed on the pairwise group dispersions. This is combined with a
+ permutation test based on the t statistic calculated on pariwise group
+ dispersions. An alternative to the classical comparison of group
+ dispersions, is to calculate Tukey's Honest Significant Differences
+ between groups, via \code{\link{TukeyHSD.betadisper}}.
+}
+\value{
+ \code{permDisper} returns a list of class \code{"permDisper"} with the
+ following components:
+
+ \item{tab}{the ANOVA table which is an object inheriting from class
+ \code{"data.frame"}.}
+ \item{pairwise}{a list with components \code{observed} and
+ \code{permuted} containing the observed and permuted p-values for
+ pairwise comparisons of group mean distances (dispersions or variances).}
+ \item{groups}{character; the levels of the grouping factor.}
+ \item{control}{a list, the result of a call to
+ \code{\link{permControl}}.}
+}
+\references{
+ Anderson, M.J. (2006) Distance-based tests for homogeneity of
+ multivariate dispersions. \emph{Biometrics} \strong{62(1)}, 245--253.
+
+ Anderson, M.J., Ellingsen, K.E. & McArdle, B.H. (2006) Multivariate
+ dispersion as a measure of beta diversity. \emph{Ecology Letters}
+ \strong{9(6)}, 683--693.
+}
+\author{Gavin L. Simpson}
+\seealso{For the main fitting function see \code{\link{betadisper}}. For
+ an alternative approach to determining which groups are more variable,
+ see \code{\link{TukeyHSD.betadisper}}.}
+\examples{
+data(varespec)
+
+## Bray-Curtis distances between samples
+dis <- vegdist(varespec)
+
+## First 16 sites grazed, remaining 8 sites ungrazed
+groups <- factor(c(rep(1,16), rep(2,8)), labels = c("grazed","ungrazed"))
+
+## Calculate multivariate dispersions
+mod <- betadisper(dis, groups)
+mod
+
+## Perform test
+anova(mod)
+
+## Permutation test for F
+permDisper(mod, pairwise = TRUE)
+
+## Tukey's Honest Significant Differences
+(mod.HSD <- TukeyHSD(mod))
+plot(mod.HSD)
+}
+\keyword{methods}
+\keyword{multivariate}
More information about the Vegan-commits
mailing list