[Vegan-commits] r286 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 25 16:58:32 CET 2008
Author: gsimpson
Date: 2008-03-25 16:58:32 +0100 (Tue, 25 Mar 2008)
New Revision: 286
Added:
pkg/R/permutest.betadisper.R
pkg/R/print.permutest.betadisper.R
pkg/man/permutest.betadisper.Rd
Removed:
pkg/R/permDisper.R
pkg/R/print.permDisper.R
pkg/man/permDisper.Rd
Log:
Changing file names due to renaming permDisper to permutest.betadisper in r285
Deleted: pkg/R/permDisper.R
===================================================================
--- pkg/R/permDisper.R 2008-03-25 15:56:10 UTC (rev 285)
+++ pkg/R/permDisper.R 2008-03-25 15:58:32 UTC (rev 286)
@@ -1,76 +0,0 @@
-`permutest.betadisper` <- function(x, 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(x, "betadisper"))
- stop("Only for class \"betadisper\"")
- nobs <- length(x$distances)
- mod <- lm(x$distances ~ x$group)
- mod.Q <- mod$qr
- p <- mod.Q$rank
- resids <- qr.resid(mod.Q, x$distances)
- res <- numeric(length = control$nperm + 1)
- res[1] <- summary(mod)$fstatistic[1]
- ## pairwise comparisons
- if(pairwise) {
- ## unique pairings
- combin <- combn(levels(x$group), 2)
- n.pairs <- ncol(combin)
- t.stats <- matrix(0, ncol = n.pairs, nrow = control$nperm + 1)
- t.stats[1,] <- apply(combn(levels(x$group), 2), 2, function(z) {
- t.statistic(x$distances[x$group == z[1]],
- x$distances[x$group == z[2]])})
- }
- for(i in seq(along = res[-1])) {
- 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)
- rss <- sum(r^2)
- 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 <- x$distance[perm][x$group == combin[1, j]]
- grp2 <- x$distance[perm][x$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(z) {
- length(x$distances[x$group == z[1]]) +
- length(x$distance[x$group == z[2]]) - 2})
- pairwise <- list(observed = 2 * pt(-abs(t.stats[1,]), df),
- permuted = apply(t.stats, 2,
- function(z) sum(abs(z) >= z[1])/length(z)))
- names(pairwise$observed) <- names(pairwise$permuted) <-
- apply(combin, 2, paste, collapse = "-")
- } else {
- pairwise <- NULL
- }
- mod.aov <- anova(x)
- 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, pairwise = pairwise,
- groups = levels(x$group), control = control)
- class(retval) <- "permutest.betadisper"
- retval
-}
Copied: pkg/R/permutest.betadisper.R (from rev 285, pkg/R/permDisper.R)
===================================================================
--- pkg/R/permutest.betadisper.R (rev 0)
+++ pkg/R/permutest.betadisper.R 2008-03-25 15:58:32 UTC (rev 286)
@@ -0,0 +1,76 @@
+`permutest.betadisper` <- function(x, 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(x, "betadisper"))
+ stop("Only for class \"betadisper\"")
+ nobs <- length(x$distances)
+ mod <- lm(x$distances ~ x$group)
+ mod.Q <- mod$qr
+ p <- mod.Q$rank
+ resids <- qr.resid(mod.Q, x$distances)
+ res <- numeric(length = control$nperm + 1)
+ res[1] <- summary(mod)$fstatistic[1]
+ ## pairwise comparisons
+ if(pairwise) {
+ ## unique pairings
+ combin <- combn(levels(x$group), 2)
+ n.pairs <- ncol(combin)
+ t.stats <- matrix(0, ncol = n.pairs, nrow = control$nperm + 1)
+ t.stats[1,] <- apply(combn(levels(x$group), 2), 2, function(z) {
+ t.statistic(x$distances[x$group == z[1]],
+ x$distances[x$group == z[2]])})
+ }
+ for(i in seq(along = res[-1])) {
+ 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)
+ rss <- sum(r^2)
+ 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 <- x$distance[perm][x$group == combin[1, j]]
+ grp2 <- x$distance[perm][x$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(z) {
+ length(x$distances[x$group == z[1]]) +
+ length(x$distance[x$group == z[2]]) - 2})
+ pairwise <- list(observed = 2 * pt(-abs(t.stats[1,]), df),
+ permuted = apply(t.stats, 2,
+ function(z) sum(abs(z) >= z[1])/length(z)))
+ names(pairwise$observed) <- names(pairwise$permuted) <-
+ apply(combin, 2, paste, collapse = "-")
+ } else {
+ pairwise <- NULL
+ }
+ mod.aov <- anova(x)
+ 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, pairwise = pairwise,
+ groups = levels(x$group), control = control)
+ class(retval) <- "permutest.betadisper"
+ retval
+}
Deleted: pkg/R/print.permDisper.R
===================================================================
--- pkg/R/print.permDisper.R 2008-03-25 15:56:10 UTC (rev 285)
+++ pkg/R/print.permDisper.R 2008-03-25 15:58:32 UTC (rev 286)
@@ -1,40 +0,0 @@
-`print.permutest.betadisper` <- function(x, digits = max(getOption("digits") - 2, 3),
- ...)
-{
- ## uses code from stats:::print.anova by R Core Development Team
- cat("\n")
- writeLines(strwrap("Permutation test for homogeneity of multivariate dispersions\n"))
- ##cat("\n")
- print(x$control)
- nc <- dim(x$tab)[2]
- cn <- colnames(x$tab)
- has.P <- substr(cn[nc], 1, 3) == "Pr("
- zap.i <- 1:(if (has.P)
- nc - 1
- else nc)
- i <- which(substr(cn, 2, 7) == " value")
- i <- c(i, which(!is.na(match(cn, "F"))))
- if (length(i))
- zap.i <- zap.i[!(zap.i %in% i)]
- tst.i <- i
- if (length(i <- grep("Df$", cn)))
- zap.i <- zap.i[!(zap.i %in% i)]
- if (length(i <- grep("N.Perm$", cn)))
- zap.i <- zap.i[!(zap.i %in% i)]
- cat("Response: Distances", sep = "\n")
- printCoefmat(x$tab, digits = digits,
- 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)
-}
Copied: pkg/R/print.permutest.betadisper.R (from rev 285, pkg/R/print.permDisper.R)
===================================================================
--- pkg/R/print.permutest.betadisper.R (rev 0)
+++ pkg/R/print.permutest.betadisper.R 2008-03-25 15:58:32 UTC (rev 286)
@@ -0,0 +1,40 @@
+`print.permutest.betadisper` <- function(x, digits = max(getOption("digits") - 2, 3),
+ ...)
+{
+ ## uses code from stats:::print.anova by R Core Development Team
+ cat("\n")
+ writeLines(strwrap("Permutation test for homogeneity of multivariate dispersions\n"))
+ ##cat("\n")
+ print(x$control)
+ nc <- dim(x$tab)[2]
+ cn <- colnames(x$tab)
+ has.P <- substr(cn[nc], 1, 3) == "Pr("
+ zap.i <- 1:(if (has.P)
+ nc - 1
+ else nc)
+ i <- which(substr(cn, 2, 7) == " value")
+ i <- c(i, which(!is.na(match(cn, "F"))))
+ if (length(i))
+ zap.i <- zap.i[!(zap.i %in% i)]
+ tst.i <- i
+ if (length(i <- grep("Df$", cn)))
+ zap.i <- zap.i[!(zap.i %in% i)]
+ if (length(i <- grep("N.Perm$", cn)))
+ zap.i <- zap.i[!(zap.i %in% i)]
+ cat("Response: Distances", sep = "\n")
+ printCoefmat(x$tab, digits = digits,
+ 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)
+}
Deleted: pkg/man/permDisper.Rd
===================================================================
--- pkg/man/permDisper.Rd 2008-03-25 15:56:10 UTC (rev 285)
+++ pkg/man/permDisper.Rd 2008-03-25 15:58:32 UTC (rev 286)
@@ -1,90 +0,0 @@
-\name{permutest.betadisper}
-\alias{permutest.betadisper}
-\alias{print.permutest.betadisper}
-\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{
-\method{permutest}{betadisper}(x, pairwise = FALSE,
- control = permControl(nperm = 999), \dots)
-}
-%- maybe also 'usage' for other objects documented here.
-\arguments{
- \item{x}{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}}}
- \item{\dots}{Arguments passed to other methods.}
-}
-\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{permutest.betadisper} 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{permutest.betadisper} returns a list of class
- \code{"permutest.betadisper"} 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
-permutest(mod, pairwise = TRUE)
-
-## Tukey's Honest Significant Differences
-(mod.HSD <- TukeyHSD(mod))
-plot(mod.HSD)
-}
-\keyword{methods}
-\keyword{multivariate}
Copied: pkg/man/permutest.betadisper.Rd (from rev 285, pkg/man/permDisper.Rd)
===================================================================
--- pkg/man/permutest.betadisper.Rd (rev 0)
+++ pkg/man/permutest.betadisper.Rd 2008-03-25 15:58:32 UTC (rev 286)
@@ -0,0 +1,90 @@
+\name{permutest.betadisper}
+\alias{permutest.betadisper}
+\alias{print.permutest.betadisper}
+\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{
+\method{permutest}{betadisper}(x, pairwise = FALSE,
+ control = permControl(nperm = 999), \dots)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{x}{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}}}
+ \item{\dots}{Arguments passed to other methods.}
+}
+\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{permutest.betadisper} 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{permutest.betadisper} returns a list of class
+ \code{"permutest.betadisper"} 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
+permutest(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