[Vegan-commits] r409 - in pkg: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 9 15:54:00 CEST 2008
Author: psolymos
Date: 2008-06-09 15:54:00 +0200 (Mon, 09 Jun 2008)
New Revision: 409
Added:
pkg/R/permat.R
pkg/R/swapcount.R
pkg/man/permat.Rd
Modified:
pkg/inst/ChangeLog
Log:
matrix permutation algorithms added
Added: pkg/R/permat.R
===================================================================
--- pkg/R/permat.R (rev 0)
+++ pkg/R/permat.R 2008-06-09 13:54:00 UTC (rev 409)
@@ -0,0 +1,177 @@
+## permatfull function
+`permatfull` <-
+function(m, fixedmar="both", reg=NULL, hab=NULL, mtype="count", times=100)
+{
+ mtype <- match.arg(mtype, c("prab", "count"))
+ count <- if (mtype == "count") TRUE else FALSE
+ fixedmar <- match.arg(fixedmar, c("none", "rows", "columns", "both"))
+ m <- as.matrix(m)
+ n.row <- nrow(m)
+ n.col <- ncol(m)
+ if (mtype == "prab") m <- matrix(as.numeric(m > 0), n.row, n.col)
+ if (is.null(reg) & is.null(hab)) str <- as.factor(rep(1, n.row))
+ if (!is.null(reg) & is.null(hab)) str <- as.factor(reg)
+ if (is.null(reg) & !is.null(hab)) str <- as.factor(hab)
+ if (!is.null(reg) & !is.null(hab)) str <- interaction(reg, hab, drop=TRUE)
+ levels(str) <- 1:length(unique(str))
+ str <- as.numeric(str)
+ nstr <- length(unique(str))
+ if (any(tapply(str,list(str),length) == 1))
+ stop("strata should contain at least 2 observations")
+ perm <- list()
+ for (k in 1:times)
+ perm[[k]] <- matrix(0, n.row, n.col)
+ for (j in 1:nstr) {
+ id <- which(str == j)
+ if (fixedmar == "none")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- matrix(sample(m[id,]), length(id), n.col)
+ else perm[[i]][id,] <- commsimulator(m[id,], method="r00")
+ if (fixedmar == "rows")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- apply(m[id,], 2, sample)
+ else perm[[i]][id,] <- commsimulator(m[id,], method="r0")
+ if (fixedmar == "cols")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- t(apply(m[id,], 1, sample))
+ else perm[[i]][id,] <- commsimulator(m[id,], method="c0")
+ if (fixedmar == "both")
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- r2dtable(1, apply(m[id,], 1, sum), apply(m[id,], 2, sum))[[1]]
+ else perm[[i]][id,] <- commsimulator(m[id,], method="quasiswap")
+ }
+ specs <- list(reg=reg, hab=hab)
+ out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
+ attr(out, "mtype") <- mtype
+ attr(out, "ptype") <- "full"
+ attr(out, "fixedmar") <- fixedmar
+ attr(out, "times") <- times
+ class(out) <- c("permat", "list")
+ return(out)
+}
+
+## permatswap function
+`permatswap` <-
+function(m, reg=NULL, hab=NULL, mtype="count", method="swap", times=100, burnin = 10000, thin = 1000)
+{
+ mtype <- match.arg(mtype, c("prab", "count"))
+ count <- if (mtype == "count") TRUE else FALSE
+ if (count) {
+ method <- match.arg(method, "swap")
+ } else {method <- match.arg(method, c("swap", "tswap", "backtrack"))}
+
+ m <- as.matrix(m)
+ n.row <- nrow(m)
+ n.col <- ncol(m)
+ if (mtype == "prab") m <- matrix(as.numeric(m > 0), n.row, n.col)
+ if (is.null(reg) & is.null(hab)) str <- as.factor(rep(1, n.row))
+ if (!is.null(reg) & is.null(hab)) str <- as.factor(reg)
+ if (is.null(reg) & !is.null(hab)) str <- as.factor(hab)
+ if (!is.null(reg) & !is.null(hab)) str <- interaction(reg, hab, drop=TRUE)
+ levels(str) <- 1:length(unique(str))
+ str <- as.numeric(str)
+ nstr <- length(unique(str))
+ if (any(tapply(str,list(str),length) == 1))
+ stop("strata should contain at least 2 observations")
+ perm <- list()
+ for (i in 1:times)
+ perm[[i]] <- matrix(0, n.row, n.col)
+
+ for (j in 1:nstr) {
+ id <- which(str == j)
+ temp <- m[id,]
+ if (count) for (k in 1:burnin) temp <- swapcount(temp)
+ else for (k in 1:burnin) temp <- commsimulator(temp, method=method)
+ for (i in 1:times)
+ if (count) perm[[i]][id,] <- swapcount(temp)
+ else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
+ }
+ specs <- list(reg=reg, hab=hab, burnin=burnin, thin=thin)
+ out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
+ attr(out, "mtype") <- mtype
+ attr(out, "ptype") <- "swap"
+ attr(out, "fixedmar") <- "both"
+ attr(out, "times") <- times
+ class(out) <- c("permat", "list")
+ return(out)
+}
+
+## S3 plot method for permat
+`plot.permat` <-
+function(x, ...)
+{
+ n <- attr(x, "times")
+ bray <- numeric(n)
+ for (i in 1:n) bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
+ plot(bray,type="n",ylim=c(0,1),ylab="Bray-Curtis dissimilarity",xlab="Runs", ...)
+ for (i in c(0.4, 0.6)) abline(i,0, lty=2, col="grey")
+ lines(bray,col="red")
+ lines(lowess(bray),col="blue",lty=2)
+ title(sub=paste("(mean = ", substitute(z, list(z=round(mean(bray),3))),
+ ", min = ", substitute(z, list(z=round(min(bray),3))),
+ ", max = ", substitute(z, list(z=round(max(bray),3))), ")", sep=""))
+ invisible(NULL)
+}
+
+## S3 print method for permat
+`print.permat` <-
+function(x, digits=3, ...)
+{
+ if (attr(x, "ptype") != "sar" & !is.null(x$specs$reg) | !is.null(x$specs$hab))
+ restr <- TRUE else restr <- FALSE
+ cat("Object of class 'permat'\n\nCall: ")
+ print(x$call)
+ cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
+ cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
+ cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
+ cat("\nSum of original matrix:", sum(x$orig))
+ cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
+ cat("\nNumber of permuted matrices:", length(x$perm),"\n")
+}
+
+## S3 summary method for permat
+`summary.permat` <-
+function(x, digits=2, ...)
+{
+ n <- attr(x, "times")
+ if (attr(x, "ptype") != "sar" & !is.null(x$specs$reg) | !is.null(x$specs$hab))
+ restr <- TRUE else restr <- FALSE
+ if (restr) {
+ if (!is.null(x$specs$reg) & is.null(x$specs$hab)) int <- x$specs$reg
+ if (is.null(x$specs$reg) & !is.null(x$specs$hab)) int <- x$specs$hab
+ if (!is.null(x$specs$reg) & !is.null(x$specs$hab))
+ int <- interaction(x$specs$reg, x$specs$hab, drop=TRUE)
+ ssum <- numeric(n)}
+ bray <- psum <- pfill <- vrow <- vcol <- numeric(n)
+ for (i in 1:n) {
+ bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
+ psum[i] <- sum(x$orig) == sum(x$perm[[i]])
+ pfill[i] <- sum(x$orig > 0) == sum(x$perm[[i]] > 0)
+ vrow[i] <- identical(rowSums(x$orig), rowSums(x$perm[[i]]))
+ vcol[i] <- identical(colSums(x$orig), colSums(x$perm[[i]]))
+ if (restr) ssum[i] <- identical(rowSums(aggregate(x$orig,list(x$specs$reg),sum)),
+ rowSums(aggregate(x$perm[[i]],list(x$specs$reg),sum)))
+ }
+ strsum <- if (restr) sum(ssum)/n else NA
+ outv <- c(sum=sum(psum)/n, fill=sum(pfill)/n, rowsums=sum(vrow)/n, colsums=sum(vcol)/n, strsum=strsum)
+
+ cat("Summary of object of class 'permat'\n\nCall: ")
+ print(x$call)
+ cat("Matrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
+ cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
+ cat("\n\nMatrix dimensions:", nrow(x$orig), "rows,", ncol(x$orig), "columns")
+ cat("\nSum of original matrix:", sum(x$orig))
+ cat("\nFill of original matrix:", round(sum(x$orig>0)/(nrow(x$orig)*ncol(x$orig)),digits))
+ cat("\nNumber of permuted matrices:", length(x$perm),"\n")
+
+ cat("\nMatrix sums retained:", round(100*outv[1], digits), "%")
+ cat("\nMatrix fill retained:", round(100*outv[2], digits), "%")
+ cat("\nRow sums retained: ", round(100*outv[3], digits), "%")
+ cat("\nColumn sums retained:", round(100*outv[4], digits), "%")
+ if (restr) cat("\nSums within strata retained:", round(100*outv[5], digits), "%")
+ cat("\n\nBray-Curtis dissimilarities among original and permuted matrices:\n")
+ print(summary(bray))
+ out <- list(bray=bray, restr=outv)
+ class(out) <- c("summary.permat", "list")
+ invisible(out)
+}
Property changes on: pkg/R/permat.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/swapcount.R
===================================================================
--- pkg/R/swapcount.R (rev 0)
+++ pkg/R/swapcount.R 2008-06-09 13:54:00 UTC (rev 409)
@@ -0,0 +1,39 @@
+`swapcount` <-
+function(m, thin = 1)
+{
+## internal, is the 2x2 matrix diagonal or anti-diagonal
+isDiag <- function(x) {
+ x<- as.vector(x)
+ x<- as.vector(x)
+ X <- as.numeric(x>0)
+ sX <- sum(X)
+ choose <- c(min(x[c(2,3)]), min(x[c(1,4)]))
+ if (sX == 4) {
+ ch <- sample(c(1,2), 1)
+ d <- choose[ch]
+ if (ch == 2) ch <- -1
+ return(d * ch)}
+ if (identical(X, c(0,1,1,0)) | identical(X, c(0,1,1,1)) | identical(X, c(1,1,1,0)))
+ return(choose[1])
+ if (identical(X, c(1,0,0,1)) | identical(X, c(1,0,1,1)) | identical(X, c(1,1,0,1)))
+ return(-choose[2])
+ if (sX < 2 | identical(X, c(0,0,1,1)) | identical(X, c(1,1,0,0)) |
+ identical(X, c(0,1,0,1)) | identical(X, c(1,0,1,0)))
+ return(0)
+ }
+ x <- as.matrix(m)
+ n.col <- ncol(x)
+ n.row <- nrow(x)
+ changed <- 0
+ while(changed < thin) {
+ ran.row <- sample(n.row, 2)
+ ran.col <- sample(n.col, 2)
+ ev <- isDiag(x[ran.row, ran.col])
+ if (ev != 0) {
+ if (identical(sum(x[ran.row, ran.col] > 0),
+ sum(x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2) > 0)))
+ x[ran.row, ran.col] <- x[ran.row, ran.col] + matrix(c(ev,-ev,-ev,ev), 2, 2)
+ changed <- changed + 1}
+ }
+ return(x)
+}
Property changes on: pkg/R/swapcount.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-06-09 08:04:12 UTC (rev 408)
+++ pkg/inst/ChangeLog 2008-06-09 13:54:00 UTC (rev 409)
@@ -4,6 +4,12 @@
Version 1.14-3 (opened June 5, 2008)
+ * permatfull, permatswap and swapcount: functions to generate
+ unrestricted and restricted null model community data matrices
+ under diferent constraints (preserving row/columnsums,
+ or incidence pattern), with print, plot and summary methods.
+ (submitted by PS)
+
* rda.default, cca.default: vegan naively used only the rank of
the ordination, but indeed, there are three cases of ranks: rank
of the ordination or number of axes (returned as rank like
Added: pkg/man/permat.Rd
===================================================================
--- pkg/man/permat.Rd (rev 0)
+++ pkg/man/permat.Rd 2008-06-09 13:54:00 UTC (rev 409)
@@ -0,0 +1,101 @@
+\name{permat}
+\alias{permatfull}
+\alias{permatswap}
+\alias{plot.permat}
+\alias{print.permat}
+\alias{summary.permat}
+\alias{swapcount}
+\title{Matrix Permutation Algorithms or Presence-Absence and Count Data}
+\description{
+Individual (for count data) or incidence (for presence-absence data) based null models can be generated for community level simulations. Options for preserving characteristics of the original matrix (sum, number of incidences, rows/columns sums) and restricted permutations (within strata based on spatial units, habitat classes or both) are discussed in the Details section. The testing of different hypothesis is separated from the null model generation, thus several tests can be applied on the same set of random matrices.
+}
+\usage{
+permatfull(m, fixedmar = "both", reg = NULL, hab = NULL, mtype = "count", times = 100)
+permatswap(m, reg = NULL, hab = NULL, mtype = "count", method = "swap", times = 100, burnin = 10000, thin = 1000)
+swapcount(m, thin = 1)
+\method{plot}{permat}(x, ...)
+\method{print}{permat}(x, digits = 3, ...)
+\method{summary}{permat}(x, digits = 2, ...)
+}
+\arguments{
+ \item{m}{a community data matrix with plots (samples) as rows and species (taxa) as columns.}
+ \item{fixedmar}{character, stating which of the row/column sums should be preserved (\code{"none", "rows", "columns", "both"}).}
+ \item{reg}{numeric vector or factor with length same as \code{nrow(m)} for grouping rows within strata (regions) for restricted permutations. Unique values or levels are used.}
+ \item{hab}{numeric vector or factor with length same as \code{nrow(m)} for grouping rows within strata (habitat classes) for restricted permutations. Unique values or levels are used.}
+ \item{mtype}{matrix data type, either \code{"count"} for count data, or \code{"prab"} for presence-absence data.}
+ \item{times}{number of permuted matrices.}
+ \item{method}{character for method used for the swap algorithm (\code{"swap", "tswap", "backtrack"}) as described for function \code{\link{commsimulator}}. If \code{mtype="count"} only \code("swap") is available using the function \code{\link{swapcount}}.}
+ \item{burnin}{number of null communities discarded before proper analysis in sequential \code{"swap", "tswap"} methods.}
+ \item{thin}{number of discarded permuted matrices between two evaluations in sequential \code{"swap", "tswap"} methods.}
+ \item{x}{object of class 'permat'}
+ \item{digits}{number of digits used for rounding.}
+ \item{\dots}{other arguments passed to methods.}
+
+}
+\details{
+Unrestricted and restricted permutations: if both \code{reg} and \code{hab} are \code{NULL}, functions perform unrestricted permutations. If either of the two is given, it is used as is for restricted permutations. If both are given, interaction is used for restricted permutations.
+
+Constraints on row/col, dill, sum
+}
+\value{
+Functions \code{permatfull} and \code{permatswap} return an object of class 'permat'.
+ \item{call}{the function call.}
+ \item{orig}{the original data matrix used for permutations.}
+ \item{perm}{a list of permuted matrices with length \code{times}.}
+ \item{specs}{a list of other specifications (variable in length, depending on the function used): \code{reg}, \code{hab}, \code{burnin}, \code{thin}.}
+
+\code{summary.permat} returns a list invisibly containing mean Bray-Curtis dissimilarities calculated pairvise among original and permuted matrices, and
+}
+\note{
+Swap methods can be very slow for large matrices.
+}
+\references{
+Original references are given on help pages of the functions used internally, listed in section 'See Also'.
+}
+\author{Peter Solymos, \email{Solymos.Peter at aotk.szie.hu}}
+\seealso{
+\code{\link{commsimulator}}, \code{\link{r2dtable}}, \code{\link{sample}}
+}
+\examples{
+library(vegan)
+
+## A simple artificial community data matrix.
+m <- matrix(c(
+ 1,3,2,0,3,1,
+ 0,2,1,0,2,1,
+ 0,0,1,2,0,3,
+ 0,0,0,1,4,3
+ ), 4, 6, byrow=TRUE)
+m
+
+## The swap algorithm for count data (1 step):
+a <- swapcount(m)
+a
+## Identiy of swapped cells:
+a != m
+
+## Using the swap algorithm to create a
+## list of permuted matrices, where
+## row/columns sums and matrix fill are preserved:
+x1 <- permatswap(m, burnin = 1000, thin = 100)
+summary(x1)
+plot(x1)
+
+## Unrestricted permutation retaining
+## row/columns sums but not matrix fill:
+x2 <- permatfull(m)
+summary(x2)
+plot(x2)
+
+## Unrestricted permutation of presence-absence type
+## retaining neither row/columns sums nor not matrix fill:
+x3 <- permatfull(m,"none", mtype="prab")
+x3$orig ## note: original matrix is binarized!
+summary(x3)
+
+## Restricted permutation,
+## check sums within strata:
+x4 <- permatfull(m,reg=c(1,1,2,2))
+summary(x4)
+}
+\keyword{multivariate}
Property changes on: pkg/man/permat.Rd
___________________________________________________________________
Name: svn:executable
+ *
More information about the Vegan-commits
mailing list