[Vegan-commits] r466 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Aug 11 23:18:46 CEST 2008
Author: psolymos
Date: 2008-08-11 23:18:46 +0200 (Mon, 11 Aug 2008)
New Revision: 466
Added:
pkg/R/adipart.R
pkg/R/permatswap.R
pkg/R/persp.tsallisaccum.R
pkg/R/plot.adipart.R
pkg/R/plot.permat.R
pkg/R/print.adipart.R
pkg/R/print.permat.R
pkg/R/print.summary.adipart.R
pkg/R/print.summary.permat.R
pkg/R/summary.adipart.R
pkg/R/summary.permat.R
pkg/R/tsallis.R
pkg/R/tsallisaccum.R
Log:
tsallis and adipart functions added to test before next major release
(adipart details undocumented)
Added: pkg/R/adipart.R
===================================================================
--- pkg/R/adipart.R (rev 0)
+++ pkg/R/adipart.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,350 @@
+## adipart function START
+adipart <-
+function(matr, strata, hclass=NULL, method="trad", index=c("richness", "shannon", "simpson"),
+ scales=seq(0, 2, 0.2), weights="unif", test=TRUE, permtype="full", times=100, crit=0.05,
+ burnin=10000, results=FALSE, ...)
+{
+ method <- match.arg(method, c("trad", "tsallis"))
+ if (method == "trad")
+ index <- match.arg(index, c("richness", "shannon", "simpson"), is.null(hclass))
+ weights <- match.arg(weights, c("unif", "prop"))
+ permtype <- match.arg(permtype, c("full", "swap"))
+ if (length(unique(strata[,1])) != 1)
+ stop("first column of strata should be uniform")
+ if (method == "tsallis" & length(scales) != 1 & !is.null(hclass))
+ stop("scales must be of length 1 if hclass is defined")
+ if (times == 0) test <- FALSE
+ if (!test & results)
+ stop("results are not produced when test is FALSE or times = 0")
+ if (!is.null(hclass) & results)
+ stop("results are not produced when hclass is defined")
+ if (inherits(matr, "matrix") | inherits(matr, "data.frame")) m <- matr
+ if (inherits(matr, "permat")) m <- matr$orig
+ if (test) {
+ if (inherits(matr, "permat")) {
+ perm <- matr$perm
+ times <- attr(matr, "times")
+ perm.done <- TRUE
+ } else {
+ perm <- NULL
+ perm.done <- FALSE}
+ } else perm <- NULL
+## internal
+nestFactor <- function(stra) {
+ nr <- nrow(stra)
+ nc <- ncol(stra)
+ for (k in 1:nc) stra[,k] <- as.numeric(stra[,k])
+ out <- stra
+ stra <- data.frame(rep(1,nrow(stra)), stra)
+ for (i in 1:nc) {
+ if (i == 1) out[,i] <- interaction(stra[,i], stra[,(i+1)], drop=TRUE)
+ else out[,i] <- interaction(out[,(i-1)], stra[,(i+1)], drop=TRUE)
+ levels(out[,i]) <- 1:nlevels(out[,i])}
+ out <- out[,c(nc:1)]
+ colnames(out) <- NULL
+ rownames(out) <- NULL
+ return(out)}
+## internal !!!f <- nestFactor(x)
+adpTrad <- function(y, f, index, weights="unif", serr=TRUE){
+ if (any(rowSums(y) == 0)) stop("empty samples not allowed")
+ ni <- length(index)
+ n <- ncol(f)
+ q <- list()
+ w <- list()
+ a <- list()
+ le <- list()
+ if (serr) sde <- list()
+ for (i in 1:n){
+ tab <- aggregate(y, by=list(f[,i]), sum)[,-1]
+ rndq <- character(3)
+ if ("richness" %in% index) {
+ Richness <- apply(tab > 0, 1, sum)
+ rndq[1] <- "Richness"
+ } else {Richness <- NULL
+ rndq[1] <- "X"}
+ if ("shannon" %in% index) {
+ Shannon <- diversity(tab, "shannon")
+ rndq[2] <- "Shannon"
+ } else {Shannon <- NULL
+ rndq[2] <- "X"}
+ if ("simpson" %in% index) {
+ Simpson <- diversity(tab, "simpson")
+ rndq[3] <- "Simpson"
+ } else {Simpson <- NULL
+ rndq[3] <- "X"}
+ rndq <- rndq[which(rndq != "X")]
+ dq.list <- list(Richness, Shannon, Simpson)
+ dq <- dqsd <- data.frame(matrix(unlist(dq.list), nrow(tab), length(rndq)))
+ colnames(dq) <- rndq
+# dq[dq == Inf | dq == -Inf] <- 0
+ le[[i]] <- length(unique(f[,i]))
+ if (weights=="prop")
+ q[[i]] <- tapply(apply(y, 1, sum), list(f[,i]), sum) / sum(y)
+ if (weights=="unif")
+ q[[i]] <- rep(1 / length(unique(f[,i])), length(unique(f[,i])))
+ if (ncol(as.matrix(dq)) != 1) {
+ for (j in 1:ni) dq[,j] <- as.matrix(dq)[,j] * q[[i]]
+ } else {dq <- dq * q[[i]]}
+ w[[i]] <- dq
+ if (ni == 1) {
+ a[[i]] <- sum(w[[i]])
+ if (serr) sde[[i]] <- sd(dqsd)
+ } else {
+ if (ncol(as.matrix(dq)) != 1) {
+ a[[i]] <- apply(as.matrix(w[[i]]), 2, sum)
+ if (serr) sde[[i]] <- apply(as.matrix(dqsd), 2, sd)
+ } else {
+ a[[i]] <- w[[i]]
+ if (serr) sde[[i]] <- sd(dqsd)}}
+ }
+ alpha <- matrix(unlist(a), n, ni, byrow=TRUE)
+ beta <- matrix(alpha[1:(n-1),], (n-1), ni, byrow=TRUE)
+ for (i in 1:(n-1))
+ beta[i, 1:ni] <- alpha[(i+1), 1:ni] - alpha[i, 1:ni]
+ colnames(alpha) <- colnames(beta) <- rndq
+ if (serr) {
+ sdev <- matrix(unlist(sde), n, ni, byrow=TRUE)
+ sqrn <- matrix(unlist(le), n, ni, byrow=FALSE)
+ sem <- sdev / sqrt(sqrn)
+# sem[is.na(sem)] <- 0
+ colnames(sem) <- rndq}
+ if (serr) {out <- list(alpha=alpha, beta=beta, sem=sem)
+ } else out <- list(alpha=alpha, beta=beta)
+ return(out)}
+## internal !!!f <- nestFactor(x)
+adpTsallis <- function(y, f, weights="unif", serr=TRUE, scales=seq(0, 2, 0.2)){
+ if (any(rowSums(y) == 0))
+ stop("empty samples not allowed")
+ n <- ncol(f)
+ ni <- length(scales)
+ q <- list()
+ w <- list()
+ a <- list()
+ le <- list()
+ if (serr) sde <- list()
+ for (i in 1:n){
+ tab <- aggregate(y, by=list(f[,i]), sum)[,-1]
+ dq <- tsallis(tab, scales=scales, norm=FALSE)
+ if (i == 1) rndq <- colnames(dq)
+ if (ni == 1) dq <- data.frame(matrix(dq, length(dq), 1))
+ dqsd <- dq
+# dq[dq == Inf | dq == -Inf] <- 0
+ le[[i]] <- length(unique(f[,i]))
+ if (weights=="prop")
+ q[[i]] <- tapply(apply(y, 1, sum), list(f[,i]), sum) / sum(y)
+ if (weights=="unif")
+ q[[i]] <- rep(1 / length(unique(f[,i])), length(unique(f[,i])))
+ if (ncol(as.matrix(dq)) != 1) {
+ for (j in 1:ni) dq[,j] <- as.matrix(dq)[,j] * q[[i]]
+ } else {dq <- dq * q[[i]]}
+ w[[i]] <- dq
+ if (any(dim(w[[i]]) != 1)) {
+ a[[i]] <- apply(w[[i]], 2, sum)
+ if (serr) sde[[i]] <- apply(dqsd, 2, sd)
+ } else {
+ if (ni == 1) {
+ a[[i]] <- apply(w[[i]], 1, sum)
+ if (serr) sde[[i]] <- sd(dqsd)}
+ if (i == n) {
+ a[[i]] <- w[[i]]
+ if (serr) sde[[i]] <- apply(as.matrix(dqsd), 1, sd)}
+ }
+ }
+ alpha <- matrix(unlist(a), n, ni, byrow=TRUE)
+ beta <- matrix(alpha[1:(n-1),], (n-1), ni, byrow=TRUE)
+ for (i in 1:(n-1))
+ beta[i, 1:ni] <- alpha[(i+1), 1:ni] - alpha[i, 1:ni]
+ colnames(alpha) <- colnames(beta) <- rndq
+ if (serr) {
+ sdev <- matrix(unlist(sde), n, ni, byrow=TRUE)
+ sqrn <- matrix(unlist(le), n, ni, byrow=FALSE)
+ sem <- sdev / sqrt(sqrn)
+# sem[is.na(sem)] <- 0
+ colnames(sem) <- rndq}
+ if (serr) {out <- list(alpha=alpha, beta=beta, sem=sem)
+ } else out <- list(alpha=alpha, beta=beta)
+ return(out)}
+## internal
+formatAdp <- function(x, style=c("alpha","beta","alpha2","beta2", "alphaH"),col.nam=NULL){
+ if (style=="alpha" | style=="alpha2" | style=="alphaH") {
+ rownames(x) <- interaction(rep("Alpha",nrow(x)), 1:nrow(x))
+ if (style=="alpha" | style=="alphaH") rownames(x)[nrow(x)] <- "Gamma"
+ if (style=="alpha2" | style=="alphaH") colnames(x) <- col.nam}
+ if (style=="beta" | style=="beta2") {
+ rownames(x) <- interaction(rep("Beta",nrow(x)), 1:nrow(x))
+ if (style=="beta2") colnames(x) <- col.nam}
+ return(x)}
+## internal
+adpPvalue <- function(obs, perm, crit=0.05, style="alpha2", col.nam=NULL){
+# perm[perm < 0] <- 0
+ Mean <- apply(perm, 1, mean)
+ obs2 <- Pval <- array(obs)
+ Sign <- (Mean - obs2) > 0
+ Abs <- abs(Mean - obs2)
+ sdp <- apply(perm, 1, sd)
+ upper <- lower <- p.val <- cl1 <- cl2 <- numeric(length(obs2))
+ adp1 <- Mean + Abs
+ adp2 <- Mean - Abs
+ for (i in 1:(nrow(obs)*ncol(obs))) {
+ upper[i] <- sum(perm[i,] >= adp1[i])
+ lower[i] <- sum(perm[i,] <= adp2[i])
+ p.val[i] <- if (Sign[i]) upper[i] else lower[i]
+ cl1[i] <- quantile(perm[i,], probs=crit/2)
+ cl2[i] <- quantile(perm[i,], probs=1-crit/2)}
+ Pval <- matrix(p.val*2 / times, nrow(obs), ncol(obs))
+# Pval[Pval > 1] <- 1
+ Mean <- matrix(Mean, nrow(obs), ncol(obs))
+ cl1 <- matrix(cl1, nrow(obs), ncol(obs))
+ cl2 <- matrix(cl2, nrow(obs), ncol(obs))
+ ses <- (obs - Mean) / matrix(sdp, nrow(obs), ncol(obs))
+ Pval <- formatAdp(Pval, style, col.nam)
+ Mean <- formatAdp(Mean, style, col.nam)
+ cl1 <- formatAdp(cl1, style, col.nam)
+ cl2 <- formatAdp(cl2, style, col.nam)
+ ses <- formatAdp(ses, style, col.nam)
+ return(list(p.value=Pval,mean=Mean,cl1=cl1,cl2=cl2,ses=ses))}
+## internal fun
+simpleCap <- function(x) {
+ s <- strsplit(x, " ")[[1]]
+ paste(toupper(substring(s, 1,1)), substring(s, 2), sep="", collapse=" ")}
+## internal
+testAdp <- function(obs, perm, f, method, results, type=c(1,2), id=NULL, burnin=NULL, ...){
+# if (type==1) id <- 1:nrow(perm[[1]])
+ pa <- matrix(NA, (nrow(obs$alpha)-1)*ncol(obs$alpha), times)
+ pb <- matrix(NA, nrow(obs$beta)*ncol(obs$alpha), times)
+ for (i in 1:times) {
+ if (perm.done) {
+ perm.i <- perm[[i]]
+ } else {
+ if (permtype == "full")
+ perm.i <- permatfull(matr, times=1, ...)$perm[[1]]
+ if (permtype == "swap") {
+ if (i == 1) {
+ perm.i <- permatswap(matr, times=1, burnin=burnin, ...)$perm[[1]]
+ } else {
+ perm.i <- permatswap(perm.i, times=1, burnin=0, ...)$perm[[1]]
+ }}}
+ if (method == "trad")
+ adp.perm <- adpTrad(perm.i[id,], f[id,], index, weights, serr=FALSE)
+ if (method == "tsallis")
+ adp.perm <- adpTsallis(perm.i[id,], f[id,], weights, serr=FALSE, scales=scales)
+ pa[,i] <- array(adp.perm$alpha[1:(nrow(adp.perm$alpha)-1),])
+ pb[,i] <- array(adp.perm$beta)}
+ obs.in.a <- obs$alpha[1:(nrow(adp.perm$alpha)-1),]
+ if (any(dim(as.matrix(obs.in.a)) == 1))
+ obs.in.a <- matrix(obs.in.a, length(obs.in.a),1)
+ obs.in.b <- obs$beta
+ if (any(dim(as.matrix(obs.in.b)) == 1))
+ obs.in.a <- matrix(obs.in.b, length(obs.in.b),1)
+ test.alpha <- adpPvalue(obs.in.a,pa,crit,"alpha2",colnames(obs$alpha))
+ test.beta <- adpPvalue(obs.in.b,pb,crit,"beta2",colnames(obs$beta))
+ if (results) res <- list(alpha=t(pa), beta=t(pb)) else res <- NULL
+ out <- list(alpha=test.alpha, beta=test.beta, res=res)
+ return(out)}
+
+## make orig
+ if (!is.null(hclass)) {
+ habnam <- unique(as.character(hclass))
+ nhab <- length(habnam)
+ habnam2 <- c(habnam, "All")
+ f2 <- data.frame(strata[,1], hclass, strata[,2:ncol(strata)])
+ obslist <- list(alpha=list(), beta=list(), sem=list())
+ for (i in 1:(nhab+1)) {
+ if (i < nhab+1) {
+ ff <- nestFactor(f2[hclass == habnam[i],])
+ mm <- m[hclass == habnam[i],]
+ } else {
+ ff <- nestFactor(f2)
+ mm <- m}
+ if (method == "trad")
+ obs <- adpTrad(mm, ff, index, weights, TRUE)
+ if (method == "tsallis")
+ obs <- adpTsallis(mm, ff, weights, TRUE, scales=scales)
+ obslist$alpha[[i]] <- obs$alpha
+ obslist$beta[[i]] <- obs$beta
+ obslist$sem[[i]] <- obs$sem
+ }
+ alpha <- formatAdp(matrix(unlist(obslist$alpha), ncol(f2), nhab+1),"alphaH",habnam2)
+ beta <- formatAdp(matrix(unlist(obslist$beta), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ sem <- formatAdp(matrix(unlist(obslist$sem), ncol(f2), nhab+1),"alphaH",habnam2)
+ obs <- list(alpha=alpha, beta=beta, sem=sem)
+ if (test) {
+ explist <- list(
+ alpha=list(p.value=list(), mean=list(), cl1=list(), cl2=list(), ses=list()),
+ beta=list(p.value=list(), mean=list(), cl1=list(), cl2=list(), ses=list()))
+ for (i in 1:(nhab+1)) {
+ if (i < nhab+1) {
+ expt <- testAdp(obs, perm, ff, method, FALSE, 2, which(hclass == habnam[i]), burnin, ...)
+ } else {
+ expt <- testAdp(obs, perm, ff, method, FALSE, 2, 1:length(hclass), burnin, ...)
+ }
+ explist$alpha$p.value[[i]] <- expt$alpha$p.value
+ explist$alpha$mean[[i]] <- expt$alpha$mean
+ explist$alpha$cl1[[i]] <- expt$alpha$cl1
+ explist$alpha$cl2[[i]] <- expt$alpha$cl2
+ explist$alpha$ses[[i]] <- expt$alpha$ses
+ explist$beta$p.value[[i]] <- expt$beta$p.value
+ explist$beta$mean[[i]] <- expt$beta$mean
+ explist$beta$cl1[[i]] <- expt$beta$cl1
+ explist$beta$cl2[[i]] <- expt$beta$cl2
+ explist$beta$ses[[i]] <- expt$beta$ses
+ }
+ alpha.p.value <- formatAdp(matrix(unlist(explist$alpha$p.value), ncol(f2)-1, nhab+1),"alpha2",habnam2)
+ alpha.mean <- formatAdp(matrix(unlist(explist$alpha$mean), ncol(f2)-1, nhab+1),"alpha2",habnam2)
+ alpha.cl1 <- formatAdp(matrix(unlist(explist$alpha$cl1), ncol(f2)-1, nhab+1),"alpha2",habnam2)
+ alpha.cl2 <- formatAdp(matrix(unlist(explist$alpha$cl2), ncol(f2)-1, nhab+1),"alpha2",habnam2)
+ alpha.ses <- formatAdp(matrix(unlist(explist$alpha$ses), ncol(f2)-1, nhab+1),"alpha2",habnam2)
+ beta.p.value <- formatAdp(matrix(unlist(explist$beta$p.value), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ beta.mean <- formatAdp(matrix(unlist(explist$beta$mean), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ beta.cl1 <- formatAdp(matrix(unlist(explist$beta$cl1), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ beta.cl2 <- formatAdp(matrix(unlist(explist$beta$cl2), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ beta.ses <- formatAdp(matrix(unlist(explist$beta$ses), ncol(f2)-1, nhab+1),"beta2",habnam2)
+ exp <- list(
+ alpha=list(p.value=alpha.p.value, mean=alpha.mean, cl1=alpha.cl1, cl2=alpha.cl2, ses=alpha.ses),
+ beta=list(p.value=beta.p.value, mean=beta.mean, cl1=beta.cl1, cl2=beta.cl2, ses=beta.ses))
+ res <- NULL
+ } else {
+ exp <- NULL
+ times <- 0
+ res <- NULL}
+ obs$alpha[(nrow(obs$alpha)-1), 1:(ncol(obs$alpha)-1)] <- NA
+ obs$beta[nrow(obs$beta), 1:(ncol(obs$beta)-1)] <- NA
+ if (test) {for (a in 1:5) {
+ exp$alpha[[a]][(nrow(obs$alpha)-1), 1:(ncol(obs$alpha)-1)] <- NA
+ exp$beta[[a]][(nrow(obs$alpha)-1), 1:(ncol(obs$alpha)-1)] <- NA
+ }}
+ index.out <- simpleCap(index)
+ design <- "twoway"
+ } else {
+ f <- nestFactor(strata)
+ if (method == "trad")
+ obs <- adpTrad(m, f, index, weights, TRUE)
+ if (method == "tsallis")
+ obs <- adpTsallis(m, f, weights, TRUE, scales=scales)
+ obs$alpha <- formatAdp(obs$alpha,"alpha")
+ obs$sem <- formatAdp(obs$sem,"alpha")
+ obs$beta <- formatAdp(obs$beta,"beta")
+ index.out <- colnames(obs$alpha)
+ design <- "oneway"
+ if (test) {
+ exp <- testAdp(obs, perm, f, method, results, 1, 1:nrow(m), burnin, ...)
+ res <- exp$res
+ exp <- exp[1:2]
+ } else {
+ exp <- NULL
+ times <- 0
+ res <- NULL}
+ }
+
+ input <- list(call=match.call(), m=m, f=nestFactor(strata), h=hclass)
+ out <- list(input=input, obs=obs, exp=exp, res=res)
+ attr(out, "method") <- method
+ attr(out, "design") <- design
+ if (method == "trad") attr(out, "index") <- index.out
+ if (method == "tsallis") attr(out, "index") <- scales
+ attr(out, "times") <- times
+ attr(out, "weights") <- weights
+ attr(out, "crit") <- crit
+ class(out) <- c("adipart", "list")
+ return(out)
+} ## adp function END
Property changes on: pkg/R/adipart.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/permatswap.R
===================================================================
--- pkg/R/permatswap.R (rev 0)
+++ pkg/R/permatswap.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,61 @@
+## permatswap function
+`permatswap` <-
+function(m, reg=NULL, hab=NULL, mtype="count", method="swap", times=100, burnin = 10000, thin = 1000)
+{
+ if (!identical(all.equal(m, round(m)), TRUE))
+ stop("function accepts only integers (counts)")
+ mtype <- match.arg(mtype, c("prab", "count"))
+ count <- mtype == "count"
+ if (count) {
+ method <- match.arg(method, "swap")
+ } else {method <- match.arg(method, c("swap", "tswap"))}
+
+ 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 <- .C("swapcount", m = as.double(temp),
+ as.integer(n.row), as.integer(n.col),
+ as.integer(1), PACKAGE = "vegan")$m
+ else
+ for (k in 1:burnin)
+ temp <- commsimulator(temp, method=method)
+ for (i in 1:times) {
+ if (count)
+ perm[[i]][id,] <- .C("swapcount",
+ m = as.double(temp),
+ as.integer(n.row),
+ as.integer(n.col),
+ as.integer(thin),
+ PACKAGE = "vegan")$m
+ else perm[[i]][id,] <- commsimulator(temp, method=method, thin=thin)
+ temp <- perm[[i]][id,]
+ }
+ }
+ 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)
+}
Property changes on: pkg/R/permatswap.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/persp.tsallisaccum.R
===================================================================
--- pkg/R/persp.tsallisaccum.R (rev 0)
+++ pkg/R/persp.tsallisaccum.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,5 @@
+persp.tsallisaccum <-
+function(x, theta = 220, phi = 15, col = heat.colors(100), zlim, ...)
+{
+persp.renyiaccum(x, theta = theta, phi = phi, col = col, zlim = zlim, ...)
+}
Property changes on: pkg/R/persp.tsallisaccum.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/plot.adipart.R
===================================================================
--- pkg/R/plot.adipart.R (rev 0)
+++ pkg/R/plot.adipart.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,147 @@
+plot.adipart <-
+function(x, rel.yax=NULL, ymax=NULL, p.legend="bottomright", ...)
+{
+# x <- object
+ if (is.null(rel.yax)) {
+ if (attr(x, "design") == "oneway") rel <- TRUE
+ if (attr(x, "design") == "twoway") rel <- FALSE
+ } else rel <- rel.yax
+ if (attr(x, "method") == "trad") {
+# main="Additive diversity partitions"
+ if (attr(x, "design") == "oneway")
+ xlab <- "Traditional diversity indices"
+ else xlab <- "Habitat classes"
+ lty <- 1}
+ if (attr(x, "method") == "tsallis") {
+# main="Additive Tsallis diversity partitions"
+ if (attr(x, "design") == "oneway")
+ xlab <- "Scale parameter"
+ else xlab <- "Habitat classes"
+ lty <- 2}
+ if (attr(x, "design") == "oneway") {
+ ylab <- if (!rel) "Diversity components" else "Relative diversity components"
+ } else {
+ if (attr(x, "method") == "tsallis") qspacer <- "q = "
+ if (attr(x, "method") == "trad") qspacer <- ""
+ spacer <- if (!rel) "D" else "Relative d"
+ ylab <- paste(spacer, "iversity partitions (", qspacer, attr(x, "index"), ")", sep="")
+ }
+ if (!is.null(x$exp)) {
+ dat <- list(obs=x$obs$alpha,
+ perm=x$exp$alpha$mean,
+ pvala=x$exp$alpha$p.value,
+ pvalb=x$exp$beta$p.value)
+ } else {
+ dat <- list(obs=x$obs$alpha,
+ perm=x$obs$alpha,
+ pvala=x$obs$alpha,
+ pvalb=x$obs$alpha)
+ dat$pvala[dat$pvala != 1] <- 1
+ dat$pvalb[dat$pvalb != 1] <- 1}
+
+ dat$perm <- t(data.frame(t(dat$perm), dat$obs[nrow(dat$obs),]))
+ dat$pvala <- t(data.frame(t(dat$pvala), rep(1, ncol(dat$obs))))
+ dat$pvalb <- t(data.frame(t(dat$pvalb), rep(1, ncol(dat$obs))))
+
+ if (attr(x, "design") == "twoway") {
+ dims <- dim(x$obs$beta)
+
+ dat$obs <- dat$obs[,c(ncol(dat$obs),1:(ncol(dat$obs)-1))]
+ dat$perm <- dat$perm[,c(ncol(dat$perm),1:(ncol(dat$perm)-1))]
+ dat$pvala <- dat$pvala[,c(ncol(dat$pvala),1:(ncol(dat$pvala)-1))]
+ dat$pvalb <- dat$pvalb[,c(ncol(dat$pvalb),1:(ncol(dat$pvalb)-1))]
+ dat$obs[dims[1], 2:dims[2]] <- dat$obs[(dims[1] + 1), 2:dims[2]]
+ dat$perm[dims[1], 2:dims[2]] <- dat$perm[(dims[1] + 1), 2:dims[2]]
+ dat$pvala[dims[1], 2:dims[2]] <- 1
+ dat$pvalb[dims[1], 2:dims[2]] <- 1}
+
+ m <- t(dat$obs)
+ mp <- t(dat$perm)
+ pdot <- t(dat$pvala < 0.05)
+ plin <- t(dat$pvalb < 0.05)
+ n <- ncol(m)
+ lwd <- 5
+ col <- "grey"
+ col2 <- "black"
+ bg <- "white"
+ pch <- 21
+ nnn <- rownames(m)
+ if (length(nnn) == 1) nnn <- c(nnn, "x")
+ fact <- factor(nnn, levels = nnn)
+ ablabs <- character(n)
+ ablabs[1] <- "alpha[1]"
+ for (i in 1:(n-1)) ablabs[(i+1)] <- paste("beta[", i, "]", sep="")
+ if (rel) for (i in 1:nrow(m)) {
+ m[i,] <- m[i,] / m[i,ncol(m)]
+ mp[i,] <- mp[i,] / mp[i,ncol(mp)]}
+ hh <- c(m[1, 1], m[1,])
+ hhh <- c(mp[1, 1], mp[1,])
+ mh <- max(max(apply(m,1,max), apply(mp,1,max))*1.05, ymax)
+ plot(fact, rep(-1000,length(fact)), type="n",
+ ylim=c(0,mh), xlim=c(-1, nrow(m)),
+# main=main,
+ xlab=xlab, ylab=ylab, ...)
+ if (length(rownames(m)) == 1)
+ fact <- factor(rownames(m), levels = rownames(m))
+ for (i in 1:n) {
+ for (j in 1:nrow(m)){
+ if (!is.null(x$exp)) {
+ if (plin[j,i]) {
+ lty <- 1
+ lwd <- 5
+ col <- "black"
+ } else {
+ lty <- 1
+ lwd <- 5
+ col <- "grey"
+ }}
+ if (i != n) arrows(j, m[j,i], j, m[j,(i+1)], code=0, lty=lty, lwd=lwd, col=col)
+ if (j == 1) {
+ arrows(-0.75, m[j,i], -1, m[j,i], code=0)
+ if (!is.null(x$exp))
+ arrows(0.25, mp[j,i], 0, mp[j,i], code=0)
+ legposo <- mean(hh[c(i, i+1)])
+ legpose <- mean(hhh[c(i, i+1)])
+ if (i == 1) {
+ text(-0.5, m[j,i], expression(alpha[1]))
+ text(-0.5, mh, "Obs")
+ if (!is.null(x$exp)) text(0.5, mp[j,i], expression(alpha[1]))
+ } else {
+ text(-0.5, legposo, substitute(beta[z], list(z=i-1)))
+ if (!is.null(x$exp)) text(0.5, legpose, substitute(beta[z], list(z=i-1)))
+ }
+ }
+ }
+ if (attr(x, "design") == "twoway") {
+ if (i != n) lines(fact, m[,i],type="l")
+ if (i != n) lines(fact, mp[,i],type="l",lty=2)
+ } else {
+ lines(fact, m[,i],type="l")
+ lines(fact, mp[,i],type="l",lty=2)}
+ }
+ for (i in n:1) {
+ for (j in 1:nrow(m)){
+ if (!is.null(x$exp)) {
+ if (pdot[j,i]) {
+ pch <- 21
+ col2 <- "white"
+ bg <- "black"
+ } else {
+ pch <- 21
+ col2 <- "black"
+ bg <- "white"}
+ }
+ points(fact[j], m[j,i],type="p", pch=pch, bg=bg, col=col2, cex=1.2)
+ }}
+
+ if (!is.null(x$exp)){
+ ptext <- paste("p <", attr(x, "crit"))
+ text(0.5, mh, "Exp")
+ lll <- strwidth(ptext)
+ legend(p.legend, lty=1, lwd=lwd-2, col=c("black","grey"),
+ legend=c(ptext, "NS"), xjust=1, text.width=lll, bty="n")
+ legend(p.legend, pch=c(19, 19), col=c("black","white"),
+ legend=c("", ""),bty="n", xjust=1, text.width=lll*1.1)
+ legend(p.legend, pch=c(21, 21), col=c("white","black"),
+ legend=c("", ""),bty="n", xjust=1, text.width=lll*1.1)}
+}
Property changes on: pkg/R/plot.adipart.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/plot.permat.R
===================================================================
--- pkg/R/plot.permat.R (rev 0)
+++ pkg/R/plot.permat.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,15 @@
+## 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",ylab="Bray-Curtis dissimilarity",xlab="Runs", ...)
+ 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)
+}
Property changes on: pkg/R/plot.permat.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/print.adipart.R
===================================================================
--- pkg/R/print.adipart.R (rev 0)
+++ pkg/R/print.adipart.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,9 @@
+print.adipart <-
+function(x, ...)
+{
+ cat("Object of class 'adipart' for additive diversity partitioning\n\nCall: ")
+ print(x$input$call)
+ cat("Method:", attr(x, "method"), "\nDesign:", attr(x, "design"), "\nWeights:", attr(x, "weights"))
+ cat("\nIndex:", attr(x, "index"))
+ cat("\nNumber of levels:", ncol(x$input$f), "\nNumber of permutations:", attr(x, "times"), "\n")
+}
Property changes on: pkg/R/print.adipart.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/print.permat.R
===================================================================
--- pkg/R/print.permat.R (rev 0)
+++ pkg/R/print.permat.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,15 @@
+## 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:", attr(x, "times"),"\n")
+}
Property changes on: pkg/R/print.permat.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/print.summary.adipart.R
===================================================================
--- pkg/R/print.summary.adipart.R (rev 0)
+++ pkg/R/print.summary.adipart.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,10 @@
+print.summary.adipart <-
+function(x, ...)
+{
+ cat("Additive diversity partitioning summary\n\nCall: ")
+ print(x$call)
+ cat("\n")
+ print(x$divres, quote=FALSE)
+ cat("---\nSignif. codes: 0 \"***\" 0.001 \"**\" 0.01 \"*\" 0.05 \".\" 0.1 \"ns\" 1\n")
+ cat("Based on", x$times, "permutations.\n")
+}
Property changes on: pkg/R/print.summary.adipart.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/print.summary.permat.R
===================================================================
--- pkg/R/print.summary.permat.R (rev 0)
+++ pkg/R/print.summary.permat.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,25 @@
+## S3 print method for summary.permat
+`print.summary.permat` <-
+function(x, digits=2, ...)
+{
+ bray <- x$bray
+ restr <- x$restr
+ test <- x$test
+ x <- x$x
+ 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:", 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), "%")
+ cat("\n\nBray-Curtis dissimilarities among original and permuted matrices:\n")
+ print(summary(bray))
+invisible(NULL)
+}
Property changes on: pkg/R/print.summary.permat.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/summary.adipart.R
===================================================================
--- pkg/R/summary.adipart.R (rev 0)
+++ pkg/R/summary.adipart.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,43 @@
+summary.adipart <-
+function(object, digits=3, ...)
+{
+## internal
+p2a <- function(x, y, digits){
+ if (length(x) != length(y)) stop("legths differ")
+ x <- round(x, digits=digits)
+ out <- y2 <- rep(NA, length(x))
+ y2[y >= 0.1] <- "ns"
+ y2[y < 0.1] <- "."
+ y2[y < 0.05] <- "*"
+ y2[y < 0.01] <- "**"
+ y2[y < 0.001] <- "***"
+ for (i in 1:length(x)) out[i] <- paste(as.character(x[i]),y2[i],sep=" ")
+ out[is.na(x)] <- ""
+ names(out) <- names(x)
+ return(out)
+}
+ x<- object
+ test <- attr(x, "times") !=0
+ obs.b <- out.b <- round(x$obs$beta[nrow(x$obs$beta):1,], digits)
+ obs.a <- out.a <- round(x$obs$alpha[(nrow(x$obs$alpha)-1):1,], digits)
+ n <- ncol(out.a)
+ Gamma <- round(x$obs$alpha[nrow(x$obs$alpha),],digits)
+ if (test) {
+ pv.b <- x$exp$beta$p.value[nrow(x$exp$beta$p.value):1,]
+ pv.a <- x$exp$alpha$p.value[nrow(x$exp$alpha$p.value):1,]
+ if (ncol(x$obs$alpha) != 1) {
+ for (i in 1:n) {
+ out.b[,i] <- p2a(obs.b[,i], pv.b[,i], digits)
+ out.a[,i] <- p2a(obs.a[,i], pv.a[,i], digits)}
+ } else {
+ out.b <- p2a(obs.b, pv.b, digits)
+ out.a <- p2a(obs.a, pv.a, digits)}
+ }
+ out <- t(data.frame(Gamma, t(out.b), t(out.a)))
+ out[is.na(out)] <- ""
+ if (ncol(x$obs$alpha) == 1)
+ colnames(out) <- colnames(x$obs$alpha)
+ output <- list(call=x$input$call, divres=out, times=attr(x, "times"))
+ class(output) <- "summary.adipart"
+return(output)
+}
Property changes on: pkg/R/summary.adipart.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/summary.permat.R
===================================================================
--- pkg/R/summary.permat.R (rev 0)
+++ pkg/R/summary.permat.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,32 @@
+## S3 summary method for permat
+`summary.permat` <-
+function(object, ...)
+{
+ 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)
+ x$perm <- NULL
+ out <- list(x=x, bray=bray, test=test, restr=restr)
+ class(out) <- c("summary.permat", "list")
+ return(out)
+}
Property changes on: pkg/R/summary.permat.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/tsallis.R
===================================================================
--- pkg/R/tsallis.R (rev 0)
+++ pkg/R/tsallis.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,35 @@
+tsallis <-
+function (x, scales = seq(0, 2, 0.2), norm=FALSE)
+{
+ x <- as.matrix(x)
+ n <- nrow(x)
+ p <- ncol(x)
+ if (p == 1) {
+ x <- t(x)
+ n <- nrow(x)
+ p <- ncol(x)
+ }
+ x <- decostand(x, "total", 1)
+ m <- length(scales)
+ result <- array(0, dim = c(n, m))
+ dimnames(result) <- list(sites = rownames(x), scale = scales)
+ for (a in 1:m) {
+ if (scales[a] != 1 && scales[a] != 0) {
+ result[, a] <- (1-(apply(x^scales[a], 1, sum)))/(scales[a] - 1)
+ }
+ else {
+ if (scales[a] == 1) result[, a] <- diversity(x, "shannon")
+ if (scales[a] == 0) result[, a] <- apply(x > 0, 1, function(y) sum(y) - 1)
+ }
+ if (norm) {
+ ST <- apply(x > 0, 1, sum)
+ if (scales[a] == 1) result[, a] <- result[, a] / log(ST)
+ else result[, a] <- result[, a] / ((ST^(1-scales[a]) - 1) / (1 - scales[a]))
+ }
+ }
+ result <- as.data.frame(result)
+ if (any(dim(result) == 1))
+ result <- unlist(result, use.names = TRUE)
+ class(result) <- c("tsallis", "renyi", class(result))
+ result
+}
Property changes on: pkg/R/tsallis.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/tsallisaccum.R
===================================================================
--- pkg/R/tsallisaccum.R (rev 0)
+++ pkg/R/tsallisaccum.R 2008-08-11 21:18:46 UTC (rev 466)
@@ -0,0 +1,45 @@
+tsallisaccum <-
+function (x, scales = seq(0, 2, 0.2), permutations = 100, raw = FALSE, ...)
+{
+ x <- as.matrix(x)
+ n <- nrow(x)
+ p <- ncol(x)
+ if (p == 1) {
+ x <- t(x)
+ n <- nrow(x)
+ p <- ncol(x)
+ }
+ m <- length(scales)
+ result <- array(dim = c(n, m, permutations))
+ dimnames(result) <- list(pooled.sites = c(1:n), scale = scales,
+ permutation = c(1:permutations))
+ for (k in 1:permutations) {
+ result[, , k] <- as.matrix(tsallis((apply(x[sample(n),
+ ], 2, cumsum)), scales = scales, ...))
+ }
+ if (raw) {
+ if (m == 1) {
+ result <- result[, 1, ]
+ }
+ }
+ else {
+ tmp <- array(dim = c(n, m, 6))
+ for (i in 1:n) {
+ for (j in 1:m) {
+ tmp[i, j, 1] <- mean(result[i, j, 1:permutations])
+ tmp[i, j, 2] <- sd(result[i, j, 1:permutations])
+ tmp[i, j, 3] <- min(result[i, j, 1:permutations])
+ tmp[i, j, 4] <- max(result[i, j, 1:permutations])
+ tmp[i, j, 5] <- quantile(result[i, j, 1:permutations],
+ 0.025)
+ tmp[i, j, 6] <- quantile(result[i, j, 1:permutations],
+ 0.975)
+ }
+ }
+ result <- tmp
+ dimnames(result) <- list(pooled.sites = c(1:n), scale = scales,
+ c("mean", "stdev", "min", "max", "Qnt 0.025", "Qnt 0.975"))
+ }
+ class(result) <- c("tsallisaccum", "renyiaccum", class(result))
+ result
+}
Property changes on: pkg/R/tsallisaccum.R
___________________________________________________________________
Name: svn:executable
+ *
More information about the Vegan-commits
mailing list