[Vegan-commits] r654 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 12 08:10:35 CET 2009
Author: psolymos
Date: 2009-01-12 08:10:35 +0100 (Mon, 12 Jan 2009)
New Revision: 654
Modified:
pkg/vegan/R/permatswap.R
pkg/vegan/R/print.summary.permat.R
pkg/vegan/R/summary.permat.R
Log:
hybrid swap method "swsh" added to permatswap to retain row/col
incidences
Modified: pkg/vegan/R/permatswap.R
===================================================================
--- pkg/vegan/R/permatswap.R 2009-01-09 06:45:39 UTC (rev 653)
+++ pkg/vegan/R/permatswap.R 2009-01-12 07:10:35 UTC (rev 654)
@@ -1,14 +1,40 @@
## permatswap function
`permatswap` <-
-function(m, method="quasiswap", strata=NULL, mtype="count", times=99, burnin = 10000, thin = 1000)
+function(m, method="quasiswap", shuffle="both", strata=NULL, mtype="count", times=99, burnin = 10000, thin = 1000)
{
+## internal function
+indshuffle <- function(x)
+{
+ N <- length(x)
+ n <- sum(x)
+ out <- numeric(N)
+ names(out) <- 1:N
+ y <- table(sample(1:N, n, replace = TRUE))
+ out[names(out) %in% names(y)] <- y
+ names(out) <- NULL
+ return(out)
+}
+bothshuffle <- function(x, y=1)
+{
+ x[x!=0] <- indshuffle(x[x!=0] - y) + y
+ return(sample(x))
+}
if (!identical(all.equal(m, round(m)), TRUE))
stop("function accepts only integers (counts)")
mtype <- match.arg(mtype, c("prab", "count"))
+ shuffle <- match.arg(shuffle, c("samp", "both"))
count <- mtype == "count"
if (count) {
- method <- match.arg(method, c("swap", "quasiswap"))
- } else {method <- match.arg(method, c("swap", "quasiswap", "tswap", "backtracking"))}
+ method <- match.arg(method, c("swap", "quasiswap", "swsh"))
+ ## warning if swapcount is to be used
+ if (method == "swap") {
+ warning("quantitative swap method may not yield random null models, use only to study its properties")
+ isSeq <- TRUE
+ } else isSeq <- FALSE
+ } else {
+ method <- match.arg(method, c("swap", "quasiswap", "tswap", "backtracking"))
+ isSeq <- method != "quasiswap"
+ }
m <- as.matrix(m)
att <- attributes(m)
@@ -36,7 +62,7 @@
temp <- m[id,]
nn.row <- nrow(m[id,])
nn.col <- ncol(m[id,])
- if (method != "quasiswap") {
+ if (isSeq) {
if (count) {
for (k in 1:burnin)
temp <- .C("swapcount", m = as.double(temp),
@@ -57,22 +83,32 @@
temp <- perm[[i]][id,]
} # for i end
} else {
- r2tabs <- r2dtable(times, rowSums(m[id,]), colSums(m[id,]))
+ if (method != "swsh") {
+ r2tabs <- r2dtable(times, rowSums(m[id,]), colSums(m[id,]))
+ } else tempPos <- temp[temp > 0]
for (i in 1:times) {
if (count) {
- ms <- sum(m[id,] > 0)
- tmp <- r2tabs[[i]]
- ## if fills are equal, no need to do it quasiswap
- if (sum(tmp > 0) != ms) {
- tmp <- .C("rswapcount",
- m = as.double(tmp),
- as.integer(nn.row),
- as.integer(nn.col),
- as.integer(ms),
- PACKAGE="vegan")$m
+ if (method != "swsh") {
+ ms <- sum(m[id,] > 0)
+ tmp <- r2tabs[[i]]
+ ## if fills are equal, no need to restore fill
+ if (sum(tmp > 0) != ms) {
+ tmp <- .C("rswapcount",
+ m = as.double(tmp),
+ as.integer(nn.row),
+ as.integer(nn.col),
+ as.integer(ms),
+ PACKAGE="vegan")$m
+ }
perm[[i]][id,] <- matrix(tmp, nrow(perm[[i]][id,]), ncol(perm[[i]][id,]))
- } else perm[[i]][id,] <- commsimulator(temp, method=method)
- }
+ } else { # method == "swsh"
+ tmp <- commsimulator(temp, method="quasiswap")
+ if (shuffle == "samp") {
+ tmp[tmp > 0] <- sample(tempPos)
+ } else tmp[tmp > 0] <- bothshuffle(tempPos)
+ perm[[i]][id,] <- tmp
+ }
+ } else perm[[i]][id,] <- commsimulator(temp, method=method)
}
thin <- burnin <- 0
}
@@ -81,9 +117,9 @@
attr(out, "mtype") <- mtype
attr(out, "ptype") <- "swap"
attr(out, "method") <- method
- attr(out, "fixedmar") <- "both"
+ attr(out, "fixedmar") <- if (method == "swsh") "none" else "both"
attr(out, "times") <- times
- attr(out, "shuffle") <- NA
+ attr(out, "shuffle") <- if (method == "swsh") shuffle else NA
attr(out, "is.strat") <- !is.null(strata)
attr(out, "strata") <- str
attr(out, "burnin") <- burnin
Modified: pkg/vegan/R/print.summary.permat.R
===================================================================
--- pkg/vegan/R/print.summary.permat.R 2009-01-09 06:45:39 UTC (rev 653)
+++ pkg/vegan/R/print.summary.permat.R 2009-01-12 07:10:35 UTC (rev 654)
@@ -28,10 +28,15 @@
cat("\nMatrix fill retained:", round(100 * sum(x$fill) / n, digits), "%")
cat("\nRow sums retained: ", round(100 * sum(x$rowsums) / n, digits), "%")
cat("\nColumn sums retained:", round(100 * sum(x$colsums) / n, digits), "%")
+ cat("\nRow incidences retained: ", round(100 * sum(x$browsums) / n, digits), "%")
+ cat("\nColumn incidences retained:", round(100 * sum(x$bcolsums) / n, digits), "%")
if (!is.null(x$strsum))
cat("\nSums within strata retained:", round(100 * sum(x$strsum) / n, digits), "%")
cat("\n\nBray-Curtis dissimilarities among original and permuted matrices:\n")
print(summary(x$bray))
-invisible(NULL)
+ cat("\nChi-squared for original matrix: ", round(x$chisq$chisq.orig, digits),
+ " (df = ", x$chisq$df, ")\n", sep = "")
+ cat("Chi-squared values among expected and permuted matrices:\n")
+ print(summary(x$chisq$chisq.perm))
+invisible(x)
}
-
Modified: pkg/vegan/R/summary.permat.R
===================================================================
--- pkg/vegan/R/summary.permat.R 2009-01-09 06:45:39 UTC (rev 653)
+++ pkg/vegan/R/summary.permat.R 2009-01-12 07:10:35 UTC (rev 654)
@@ -8,6 +8,8 @@
fi <- sum(x$orig > 0)
rs <- rowSums(x$orig)
cs <- colSums(x$orig)
+ rb <- rowSums(x$orig > 0)
+ cb <- colSums(x$orig > 0)
nr <- nrow(x$orig)
nc <- ncol(x$orig)
bray <- sapply(x$perm, function(z) sum(abs(x$orig - z)) / sum(x$orig + z))
@@ -15,6 +17,8 @@
pfill <- sapply(x$perm, function(z) fi == sum(z > 0))
vrow <- sapply(x$perm, function(z) sum(rs == rowSums(z)) == nr)
vcol <- sapply(x$perm, function(z) sum(cs == colSums(z)) == nc)
+ brow <- sapply(x$perm, function(z) sum(rb == rowSums(z > 0)) == nr)
+ bcol <- sapply(x$perm, function(z) sum(cb == colSums(z > 0)) == nc)
if (attr(x, "is.strat")) {
int <- attr(x, "strata")
nlev <- length(unique(int))
@@ -22,8 +26,15 @@
ssum <- sapply(x$perm, function(z)
sum(rsagg == rowSums(aggregate(z, list(int), sum)[,-1])) == nlev)
} else ssum <- NULL
+ ## Chisq
+ E <- rs %o% cs / ss
+ chisq <- list(
+ chisq.orig = sum((x$orig - E)^2 / E),
+ df = (nr - 1) * (nc - 1),
+ chisq.perm = sapply(x$perm, function(z) sum((z - E)^2 / E)))
x$perm <- NULL
- out <- list(x=x, bray=bray, sum=psum, fill=pfill, rowsums=vrow, colsums=vcol, strsum=ssum)
+ out <- list(x=x, bray=bray, chisq=chisq, sum=psum, fill=pfill, rowsums=vrow, colsums=vcol,
+ browsums=brow, bcolsums=bcol, strsum=ssum)
class(out) <- c("summary.permat", "list")
return(out)
}
More information about the Vegan-commits
mailing list