[Vegan-commits] r620 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 7 20:02:09 CET 2008
Author: psolymos
Date: 2008-12-07 20:02:09 +0100 (Sun, 07 Dec 2008)
New Revision: 620
Modified:
pkg/vegan/R/permatfull.R
pkg/vegan/R/permatswap.R
pkg/vegan/R/plot.permat.R
pkg/vegan/R/print.permat.R
pkg/vegan/R/print.summary.permat.R
pkg/vegan/R/summary.permat.R
Log:
this is the big permat update
Modified: pkg/vegan/R/permatfull.R
===================================================================
--- pkg/vegan/R/permatfull.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/permatfull.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -1,6 +1,6 @@
## permatfull function
`permatfull` <-
-function(m, fixedmar="both", shuffle="ind", reg=NULL, hab=NULL, mtype="count", times=100)
+function(m, fixedmar="both", shuffle="ind", strata=NULL, mtype="count", times=99)
{
## internal function
indshuffle <- function(x)
@@ -33,10 +33,11 @@
n.row <- nrow(m)
n.col <- ncol(m)
if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
- 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)
+
+ if (is.null(strata))
+ str <- as.factor(rep(1, n.row))
+ else str <- as.factor(strata)[drop = TRUE]
+
levels(str) <- 1:length(unique(str))
str <- as.numeric(str)
nstr <- length(unique(str))
@@ -67,14 +68,17 @@
}
if (fixedmar == "both")
shuffle <- NA
- specs <- list(reg=reg, hab=hab)
- out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
+ out <- list(call=match.call(), orig=m, perm=perm)
attr(out, "mtype") <- mtype
attr(out, "ptype") <- "full"
attr(out, "method") <- NA
attr(out, "fixedmar") <- fixedmar
attr(out, "times") <- times
attr(out, "shuffle") <- shuffle
+ attr(out, "is.strat") <- !is.null(strata)
+ attr(out, "strata") <- str
+ attr(out, "burnin") <- NA
+ attr(out, "thin") <- NA
class(out) <- c("permat", "list")
return(out)
}
Modified: pkg/vegan/R/permatswap.R
===================================================================
--- pkg/vegan/R/permatswap.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/permatswap.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -1,6 +1,6 @@
## permatswap function
`permatswap` <-
-function(m, method="quasiswap", reg=NULL, hab=NULL, mtype="count", times=100, burnin = 10000, thin = 1000)
+function(m, method="quasiswap", strata=NULL, mtype="count", times=99, burnin = 10000, thin = 1000)
{
if (!identical(all.equal(m, round(m)), TRUE))
stop("function accepts only integers (counts)")
@@ -15,10 +15,11 @@
n.row <- nrow(m)
n.col <- ncol(m)
if (mtype == "prab") m <- ifelse(m > 0, 1, 0)
- 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)
+
+ if (is.null(strata))
+ str <- as.factor(rep(1, n.row))
+ else str <- as.factor(strata)[drop = TRUE]
+
levels(str) <- 1:length(unique(str))
str <- as.numeric(str)
nstr <- length(unique(str))
@@ -75,14 +76,17 @@
thin <- burnin <- 0
}
} # for j end
- specs <- list(reg=reg, hab=hab, burnin=burnin, thin=thin)
- out <- list(call=match.call(), orig=m, perm=perm, specs=specs)
+ out <- list(call=match.call(), orig=m, perm=perm)
attr(out, "mtype") <- mtype
attr(out, "ptype") <- "swap"
attr(out, "method") <- method
attr(out, "fixedmar") <- "both"
attr(out, "times") <- times
attr(out, "shuffle") <- NA
+ attr(out, "is.strat") <- !is.null(strata)
+ attr(out, "strata") <- str
+ attr(out, "burnin") <- burnin
+ attr(out, "thin") <- thin
class(out) <- c("permat", "list")
return(out)
}
Modified: pkg/vegan/R/plot.permat.R
===================================================================
--- pkg/vegan/R/plot.permat.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/plot.permat.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -4,13 +4,15 @@
{
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]])
+# for (i in 1:n) bray[i] <- sum(abs(x$orig-x$perm[[i]]))/sum(x$orig+x$perm[[i]])
+ bray <- sapply(x$perm, function(z) sum(abs(x$orig - z)) / sum(x$orig + z))
if (plot) {
plot(bray,type="n",ylab=ylab,xlab=xlab, ...)
lines(bray,col=col[1], lty=lty[1])
- lines(lowess(bray),col=col[2], lty=lty[2])}
- if (text) 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=""))
+ lines(lowess(bray),col=col[2], lty=lty[2])
+ if (text) 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(bray)
}
Modified: pkg/vegan/R/print.permat.R
===================================================================
--- pkg/vegan/R/print.permat.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/print.permat.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -2,21 +2,26 @@
`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("Object of class 'permat'\n")
cat("\nMatrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
- if (attr(x, "ptype") == "swap")
- cat("\nMethod:", attr(x, "method"))
- cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
+ if (attr(x, "ptype") == "swap") {
+ cat("\nMethod: ", attr(x, "method"), sep = "")
+ if (attr(x, "method") != "quasiswap") {
+ cat(", burnin: ", attr(x, "burnin"), sep = "")
+ cat(", thin: ", attr(x, "thin"), sep = "")
+ }
+ }
+ cat("\nRestricted:", attr(x, "is.strat"), "\nFixed margins:", attr(x, "fixedmar"))
if (!is.na(attr(x, "shuffle"))) {
if (attr(x, "shuffle")=="ind") cat("\nIndividuals")
if (attr(x, "shuffle")=="samp") cat("\nSamples")
if (attr(x, "shuffle")=="both") cat("\nIndividuals and samples")
- cat(" are shuffled")}
- 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:", attr(x, "times"),"\n")
+ cat(" are shuffled")
+ }
+ cat("\n")
+ invisible(x)
+# 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:", attr(x, "times"),"\n")
}
Modified: pkg/vegan/R/print.summary.permat.R
===================================================================
--- pkg/vegan/R/print.summary.permat.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/print.summary.permat.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -2,31 +2,36 @@
`print.summary.permat` <-
function(x, digits=2, ...)
{
- bray <- x$bray
- restr <- x$restr
- test <- x$test
- x <- x$x
+ n <- attr(x$x, "times")
cat("Summary of object of class 'permat'\n\nCall: ")
- print(x$call)
- cat("\nMatrix type:", attr(x, "mtype"), "\nPermutation type:", attr(x, "ptype"))
- if (attr(x, "ptype") == "swap")
- cat("\nMethod:", attr(x, "method"))
- cat("\nRestricted:", restr, "\nFixed margins:", attr(x, "fixedmar"))
- if (!is.na(attr(x, "shuffle"))) {
- if (attr(x, "shuffle")=="ind") cat("\nIndividuals")
- if (attr(x, "shuffle")=="samp") cat("\nSamples")
- if (attr(x, "shuffle")=="both") cat("\nIndividuals and samples")
- cat(" are shuffled")}
- 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:", attr(x, "times"),"\n")
- cat("\nMatrix sums retained:", round(100*test[1], digits), "%")
- cat("\nMatrix fill retained:", round(100*test[2], digits), "%")
- cat("\nRow sums retained: ", round(100*test[3], digits), "%")
- cat("\nColumn sums retained:", round(100*test[4], digits), "%")
- if (restr) cat("\nSums within strata retained:", round(100*test[5], digits), "%")
+ print(x$x$call)
+ cat("\nMatrix type:", attr(x$x, "mtype"), "\nPermutation type:", attr(x$x, "ptype"))
+ if (attr(x$x, "ptype") == "swap") {
+ cat("\nMethod: ", attr(x$x, "method"), sep = "")
+ if (attr(x$x, "method") != "quasiswap") {
+ cat(", burnin: ", attr(x$x, "burnin"), sep = "")
+ cat(", thin: ", attr(x$x, "thin"), sep = "")
+ }
+ }
+ cat("\nRestricted:", attr(x$x, "is.strat"), "\nFixed margins:", attr(x$x, "fixedmar"))
+ if (!is.na(attr(x$x, "shuffle"))) {
+ if (attr(x$x, "shuffle")=="ind") cat("\nIndividuals")
+ if (attr(x$x, "shuffle")=="samp") cat("\nSamples")
+ if (attr(x$x, "shuffle")=="both") cat("\nIndividuals and samples")
+ cat(" are shuffled")
+ }
+ cat("\n\nMatrix dimensions:", nrow(x$x$orig), "rows,", ncol(x$x$orig), "columns")
+ cat("\nSum of original matrix:", sum(x$x$orig))
+ cat("\nFill of original matrix:", round(sum(x$x$orig>0)/(nrow(x$x$orig)*ncol(x$x$orig)),digits))
+ cat("\nNumber of permuted matrices:", n,"\n")
+ cat("\nMatrix sums retained:", round(100 * sum(x$sum) / n, digits), "%")
+ 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), "%")
+ if (!is.null(x$restr))
+ 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(bray))
+ print(summary(x$bray))
invisible(NULL)
}
+
Modified: pkg/vegan/R/summary.permat.R
===================================================================
--- pkg/vegan/R/summary.permat.R 2008-12-07 19:01:17 UTC (rev 619)
+++ pkg/vegan/R/summary.permat.R 2008-12-07 19:02:09 UTC (rev 620)
@@ -4,29 +4,26 @@
{
x <- object
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)
- nlev <- length(unique(int))
- 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] <- sum(rowSums(x$orig) == rowSums(x$perm[[i]])) == nrow(x$orig)
- vcol[i] <- sum(colSums(x$orig) == colSums(x$perm[[i]])) == ncol(x$orig)
- if (restr) ssum[i] <- {sum(rowSums(aggregate(x$orig,list(int),sum)[,-1]) ==
- rowSums(aggregate(x$perm[[i]],list(int),sum)[,-1])) == nlev}
- }
- strsum <- if (restr) sum(ssum)/n else NA
- test <- c(sum=sum(psum)/n, fill=sum(pfill)/n, rowsums=sum(vrow)/n, colsums=sum(vcol)/n, strsum=strsum)
+ ss <- sum(x$orig)
+ fi <- sum(x$orig > 0)
+ rs <- rowSums(x$orig)
+ cs <- colSums(x$orig)
+ nr <- nrow(x$orig)
+ nc <- ncol(x$orig)
+ bray <- sapply(x$perm, function(z) sum(abs(x$orig - z)) / sum(x$orig + z))
+ psum <- sapply(x$perm, function(z) ss == sum(z))
+ 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)
+ if (attr(x, "is.strat")) {
+ int <- attr(x, "strata")
+ nlev <- length(unique(int))
+ rsagg <- rowSums(aggregate(x$orig, list(int), sum)[,-1])
+ ssum <- sapply(x$perm, function(z)
+ sum(rsagg == rowSums(aggregate(z, list(int), sum)[,-1])) == nlev)
+ } else ssum <- NULL
x$perm <- NULL
- out <- list(x=x, bray=bray, test=test, restr=restr)
+ out <- list(x=x, bray=bray, sum=psum, fill=pfill, rowsums=vrow, colsums=vcol, strsum=ssum)
class(out) <- c("summary.permat", "list")
return(out)
}
More information about the Vegan-commits
mailing list