[Vegan-commits] r828 - in pkg/vegan: R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue May 19 18:00:53 CEST 2009


Author: psolymos
Date: 2009-05-19 18:00:52 +0200 (Tue, 19 May 2009)
New Revision: 828

Modified:
   pkg/vegan/R/adipart.R
   pkg/vegan/R/multipart.R
   pkg/vegan/R/print.adipart.R
   pkg/vegan/R/print.multipart.R
   pkg/vegan/man/adipart.Rd
Log:
nsimul=0 option and relative argument added

Modified: pkg/vegan/R/adipart.R
===================================================================
--- pkg/vegan/R/adipart.R	2009-05-19 08:43:01 UTC (rev 827)
+++ pkg/vegan/R/adipart.R	2009-05-19 16:00:52 UTC (rev 828)
@@ -13,6 +13,10 @@
     nlevs <- length(tlab)
     if (nlevs < 2)
         stop("provide at least two level hierarchy")
+    if (any(rowSums(lhs) == 0))
+        stop("data matrix contains empty rows")
+    if (any(lhs < 0))
+        stop("data matrix contains negative entries")
 
     ## part check proper design of the model frame
     noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
@@ -77,8 +81,15 @@
         b <- sapply(2:nlevs, function(i) a[i] - a[(i-1)])
         c(a, b)
     }
-    sim <- oecosimu(lhs, wdivfun, method = "permat", nsimul=nsimul,
-        burnin=control$burnin, thin=control$thin, control=control)
+    if (nsimul > 0) {
+        sim <- oecosimu(lhs, wdivfun, method = "permat", nsimul=nsimul,
+            burnin=control$burnin, thin=control$thin, control=control)
+        } else {
+            sim <- wdivfun(lhs)
+            tmp <- rep(NA, length(sim))
+            sim <- list(statistic = sim,
+                oecosimu = list(z = tmp, pval = tmp, method = NA, statistic = sim))
+        }
     nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
         paste("beta", 1:(nlevs-1), sep="."))
     names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam

Modified: pkg/vegan/R/multipart.R
===================================================================
--- pkg/vegan/R/multipart.R	2009-05-19 08:43:01 UTC (rev 827)
+++ pkg/vegan/R/multipart.R	2009-05-19 16:00:52 UTC (rev 828)
@@ -1,6 +1,6 @@
 multipart <-
 function(formula, data, index=c("renyi", "tsallis"), scales = 1,
-    global = TRUE, nsimul=99, control, ...)
+    global = FALSE, relative = FALSE, nsimul=99, control, ...)
 {
     if (length(scales) > 1)
         stop("length of 'scales' must be 1")
@@ -15,6 +15,10 @@
     nlevs <- length(tlab)
     if (nlevs < 2)
         stop("provide at least two level hierarchy")
+    if (any(rowSums(lhs) == 0))
+        stop("data matrix contains empty rows")
+    if (any(lhs < 0))
+        stop("data matrix contains negative entries")
 
     ## part check proper design of the model frame
     noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
@@ -76,9 +80,12 @@
             } else {
                 tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
             }
-            a <- sapply(1:nlevs, function(i) mean(divfun(tmp[[i]]), na.rm=TRUE))
+            raw <- sapply(1:nlevs, function(i) divfun(tmp[[i]]))
+            a <- sapply(raw, mean)
             G <- a[nlevs]
             b <- sapply(1:(nlevs-1), function(i) G / a[i])
+            if (relative)
+                b <- b / sapply(raw[1:(nlevs-1)], length)
             c(a, b)
         }
     } else {
@@ -95,12 +102,22 @@
                         mean(a[[i]][id[[i]]==ii])
                     })
                 })
-            b <- sapply(1:(nlevs-1), function(i) mean(a[[(i+1)]] / am[[i]]))
-            c(sapply(a, mean), b)
+            b <- lapply(1:(nlevs-1), function(i) a[[(i+1)]] / am[[i]])
+            bmax <- lapply(id, function(i) table(i))
+            if (relative)
+                b <- lapply(1:(nlevs-1), function(i) b[[i]] / bmax[[i]])
+            c(sapply(a, mean), sapply(b, mean))
         }
     }
-    sim <- oecosimu(lhs, wdivfun, method = "permat", nsimul=nsimul,
-        burnin=control$burnin, thin=control$thin, control=control)
+    if (nsimul > 0) {
+            sim <- oecosimu(lhs, wdivfun, method = "permat", nsimul=nsimul,
+                burnin=control$burnin, thin=control$thin, control=control)
+        } else {
+            sim <- wdivfun(lhs)
+            tmp <- rep(NA, length(sim))
+            sim <- list(statistic = sim,
+                oecosimu = list(z = tmp, pval = tmp, method = NA, statistic = sim))
+        }
     nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
         paste("beta", 1:(nlevs-1), sep="."))
     names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam

Modified: pkg/vegan/R/print.adipart.R
===================================================================
--- pkg/vegan/R/print.adipart.R	2009-05-19 08:43:01 UTC (rev 827)
+++ pkg/vegan/R/print.adipart.R	2009-05-19 16:00:52 UTC (rev 828)
@@ -1,7 +1,9 @@
 print.adipart <-
 function(x, ...)
 {
-    cat("adipart with", ncol(x$oecosimu$simulated), "simulations\n")
+    n <- if (is.null(x$oecosimu$simulated))
+        0 else ncol(x$oecosimu$simulated)
+    cat("adipart with", n, "simulations\n")
     att <- attributes(x)
     att$names <- att$call <- att$class <- att$n.levels <- att$terms <- att$model <- NULL
     cat("with", paste(names(att), att, collapse=", "))
@@ -12,7 +14,12 @@
         NextMethod("print", x)
         cat("\n")
     }
-    qu <- apply(x$oecosimu$simulated, 1, quantile, probs=c(0.025, 0.5, 0.975))
+    if (!is.null(x$oecosimu$simulated)) {
+        tmp <- x$oecosimu$simulated
+    } else {
+        tmp <- data.matrix(x$oecosimu$statistic)
+    }
+    qu <- apply(tmp, 1, quantile, probs=c(0.025, 0.5, 0.975))
     m <- cbind("statistic" = x$oecosimu$statistic,
                "z" = x$oecosimu$z, t(qu),
                "Pr(sim.)"=x$oecosimu$pval)

Modified: pkg/vegan/R/print.multipart.R
===================================================================
--- pkg/vegan/R/print.multipart.R	2009-05-19 08:43:01 UTC (rev 827)
+++ pkg/vegan/R/print.multipart.R	2009-05-19 16:00:52 UTC (rev 828)
@@ -1,7 +1,9 @@
 print.multipart <-
 function(x, ...)
 {
-    cat("multipart with", ncol(x$oecosimu$simulated), "simulations\n")
+    n <- if (is.null(x$oecosimu$simulated))
+        0 else ncol(x$oecosimu$simulated)
+    cat("multipart with", n, "simulations\n")
     att <- attributes(x)
     att$names <- att$call <- att$class <- att$n.levels <- att$terms <- att$model <- NULL
     cat("with", paste(names(att), att, collapse=", "))
@@ -11,7 +13,12 @@
         NextMethod("print", x)
         cat("\n")
     }
-    qu <- apply(x$oecosimu$simulated, 1, quantile, probs=c(0.025, 0.5, 0.975))
+    if (!is.null(x$oecosimu$simulated)) {
+        tmp <- x$oecosimu$simulated
+    } else {
+        tmp <- data.matrix(x$oecosimu$statistic)
+    }
+    qu <- apply(tmp, 1, quantile, probs=c(0.025, 0.5, 0.975))
     m <- cbind("statistic" = x$oecosimu$statistic,
                "z" = x$oecosimu$z, t(qu),
                "Pr(sim.)"=x$oecosimu$pval)

Modified: pkg/vegan/man/adipart.Rd
===================================================================
--- pkg/vegan/man/adipart.Rd	2009-05-19 08:43:01 UTC (rev 827)
+++ pkg/vegan/man/adipart.Rd	2009-05-19 16:00:52 UTC (rev 828)
@@ -18,7 +18,7 @@
 adipart(formula, data, index=c("richness", "shannon", "simpson"),
     weights=c("unif", "prop"), relative = FALSE, nsimul=99, control, ...)
 multipart(formula, data, index=c("renyi", "tsallis"), scales = 1,
-    global = TRUE, nsimul=99, control, ...)
+    global = FALSE, relative = FALSE, nsimul=99, control, ...)
 hiersimu(formula, data, FUN, location = c("mean", "median"),
     relative = FALSE, drop.highest = FALSE, nsimul=99, control, ...)
 \method{print}{adipart}(x, ...)
@@ -39,11 +39,14 @@
     weighting proportional to sample abundances to use in weighted averaging of individual 
     alpha values within strata of a given level of the sampling hierarchy.}
   \item{relative}{Logical, if \code{TRUE} then alpha and beta diversity values are given 
-    relative to the value of gamma. This can be useful when comparing different indices.}
+    relative to the value of gamma for function \code{adipart}. For function \code{multipart},
+    it sets the standardization of the beta diversity value with its maximum (see Details).}
   \item{scales}{Numeroc, of length 1, the order of the generalized diversity index 
     to be used.}
   \item{global}{Logical, indicates the calculation of beta values, see Details.}
-  \item{nsimul}{Number of permutation to use if \code{matr} is not of class 'permat'.}
+  \item{nsimul}{Number of permutation to use if \code{matr} is not of class 'permat'.
+    If \code{nsimul = 0}, only the \code{FUN} argument is evaluated. It is thus possible
+    to reuse the statistic values without using a null model.}
   \item{control}{A list of arguments passed to quantitative permutation
           algorithms. If missing, the function  'permat.control' is used.}
   \item{FUN}{A function to be used by \code{hiersimu}.}
@@ -106,6 +109,10 @@
 level \eqn{i} is the mean of the beta values of the clusters at that level,
 \eqn{\beta_{i} = mean(\beta_{ij})}.
 
+If \code{relative = TRUE} for \code{multipart}, the respective beta diversity values are
+standardized by their maximum expected values (\eqn{mean(\beta_{ij}) / \beta_{max,ij}})
+given as \eqn{\beta_{max,ij} = n_{j}} (the number of lower level units in a given cluster \eqn{j}).
+
 The expected diversity components are calculated \code{nsimul} times by individual based 
 randomisation of the community data matrix. This is done by the \code{\link{permatfull}} 
 and \code{\link{permatswap}} functions, and properties of the null model can be set by 



More information about the Vegan-commits mailing list