[Vegan-commits] r954 - in branches/1.15: R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Aug 26 17:17:42 CEST 2009


Author: jarioksa
Date: 2009-08-26 17:17:42 +0200 (Wed, 26 Aug 2009)
New Revision: 954

Added:
   branches/1.15/R/estaccumR.R
   branches/1.15/R/plot.poolaccum.R
   branches/1.15/R/poolaccum.R
   branches/1.15/R/print.poolaccum.R
   branches/1.15/R/summary.poolaccum.R
Modified:
   branches/1.15/R/specaccum.R
   branches/1.15/R/specpool.R
   branches/1.15/R/specpool2vect.R
   branches/1.15/inst/ChangeLog
   branches/1.15/man/specaccum.Rd
   branches/1.15/man/specpool.Rd
Log:
merge poolaccum & estaccumR to branches/1.15 (r913 to 923)

Copied: branches/1.15/R/estaccumR.R (from rev 923, pkg/vegan/R/estaccumR.R)
===================================================================
--- branches/1.15/R/estaccumR.R	                        (rev 0)
+++ branches/1.15/R/estaccumR.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -0,0 +1,21 @@
+##" Individual based accumulation model. Similar to poolaccum but uses
+##estimateR. Inherits from "poolaccum" class and uses its methods.
+`estaccumR` <-
+    function(x, permutations = 100)
+{
+    n <- nrow(x)
+    N <- seq_len(n)
+    S <- chao <- ace <- matrix(0, nrow = n, ncol = permutations)
+    for (i in 1:permutations) {
+        take <- sample(n)
+        tmp <- estimateR(apply(x[take,], 2, cumsum))
+        S[,i] <- tmp[1,]
+        chao[,i] <- tmp[2,]
+        ace[, i] <- tmp[4,]
+    }
+    means <- cbind(N = N, S = rowMeans(S), Chao = rowMeans(chao),
+                   ACE = rowMeans(ace))
+    out <- list(S = S, chao = chao, ace = ace, N = N, means = means)
+    class(out) <- c("estaccumR", "poolaccum")
+    out
+}

Copied: branches/1.15/R/plot.poolaccum.R (from rev 923, pkg/vegan/R/plot.poolaccum.R)
===================================================================
--- branches/1.15/R/plot.poolaccum.R	                        (rev 0)
+++ branches/1.15/R/plot.poolaccum.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -0,0 +1,14 @@
+`plot.poolaccum` <-
+    function(x, alpha = 0.05, type = c("l","g"), ...)
+{
+    require(lattice) || stop("Needs package 'lattice'")
+    m <- summary(x, alpha = alpha, ...)
+    n <- nrow(m[[1]])
+    Size <- as.vector(sapply(m, function(x) c(x[,1], x[,1], rev(x[,1]))))
+    Richness <- as.vector(sapply(m, function(x) c(x[,2], x[,3], rev(x[,4]))) )
+    indnames <- as.character(sapply(m, function(x) colnames(x[,2, drop=FALSE])))
+    Index <- factor(rep(indnames, each = 3 * n), levels = indnames)
+    lintype <- rep(c(rep("aver", n), rep("envel", 2*n)), length=length(Size))
+    xyplot(Richness ~  Size | Index, as.table = TRUE, groups = lintype,
+           type = type, ...)
+}

Copied: branches/1.15/R/poolaccum.R (from rev 923, pkg/vegan/R/poolaccum.R)
===================================================================
--- branches/1.15/R/poolaccum.R	                        (rev 0)
+++ branches/1.15/R/poolaccum.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -0,0 +1,37 @@
+`poolaccum` <-
+    function(x, permutations = 100, minsize = 3)
+{
+    x <- as.matrix(x)
+    n <- nrow(x)
+    m <- ncol(x)
+    N <- seq_len(n)
+    S <- chao <- boot <- jack1 <- jack2 <-
+        matrix(0, nrow=n, ncol=permutations)
+    ## specpool() is slow, but the vectorized versions below are
+    ## pretty fast
+    for (i in 1:permutations) {
+        take <- sample(n)
+        tmp <- apply(x[take,] > 0, 2, cumsum)
+        S[,i] <- rowSums(tmp > 0)
+        ## All-zero species are taken as *known* to be missing in
+        ## subsamples, and in the following we subtract them (as
+        ## 2*S-m) from the bootstrap samples to give a more unbiased
+        ## estimate.
+        boot[,i] <- 2*S[,i] - m + rowSums(exp(sweep(log1p(-sweep(tmp, 1, N, "/")), 1, N, "*") ))
+        a1 <- rowSums(tmp == 1)
+        a2 <- rowSums(tmp == 2)
+        chao[, i] <- S[,i] + ifelse(a2 > 0, a1*a1/2/a2, 0)
+        jack1[,i] <- S[,i] + a1 * (N-1)/N
+        jack2[,i] <- S[,i] + a1*(2*N-3)/N - a2*(N-2)^2/N/(N-1)
+    }
+    means <- cbind(`N` = N, `S` = rowMeans(S), `Chao` =  rowMeans(chao),
+                   `Jackknife 1` = rowMeans(jack1),
+                   `Jackknife 2` = rowMeans(jack2),
+                   `Bootstrap` = rowMeans(boot))
+    take <- N >= minsize
+    out <- list(S = S[take,], chao = chao[take,], jack1 = jack1[take,],
+                jack2 = jack2[take,], boot = boot[take,], N = N[take],
+                means = means[take,])
+    class(out) <- "poolaccum"
+    out
+}

Copied: branches/1.15/R/print.poolaccum.R (from rev 923, pkg/vegan/R/print.poolaccum.R)
===================================================================
--- branches/1.15/R/print.poolaccum.R	                        (rev 0)
+++ branches/1.15/R/print.poolaccum.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -0,0 +1,7 @@
+`print.poolaccum` <-
+    function(x, ...)
+{
+    rownames(x$means) <- rep("", nrow(x$means))
+    print(x$means, ...)
+    invisible(x)
+}

Modified: branches/1.15/R/specaccum.R
===================================================================
--- branches/1.15/R/specaccum.R	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/R/specaccum.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -1,6 +1,6 @@
 `specaccum` <-
     function (comm, method = "exact", permutations = 100, conditioned=TRUE,
-              gamma="Jack.1", ...)
+              gamma="jack1", ...)
 {
     x <- comm
     x <- as.matrix(x)

Modified: branches/1.15/R/specpool.R
===================================================================
--- branches/1.15/R/specpool.R	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/R/specpool.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -48,9 +48,9 @@
         Zp <- (crossprod(X == 0)/n)^n - outer(pn, pn, "*")
         var.boot[is] <- sum(pn * (1 - pn)) + 2 * sum(Zp[lower.tri(Zp)])
     }
-    out <- list(Species = S, Chao = chao,  Chao.SE = sqrt(var.chao), 
-                Jack.1 = jack.1, Jack1.SE = sqrt(var.jack1), Jack.2 = jack.2, 
-                Boot = bootS, Boot.SE = sqrt(var.boot), n = as.vector(groups))
+    out <- list(Species = S, chao = chao,  chao.se = sqrt(var.chao), 
+                jack1 = jack.1, jack1.se = sqrt(var.jack1), jack2 = jack.2, 
+                boot = bootS, boot.se = sqrt(var.boot), n = as.vector(groups))
     out <- as.data.frame(out)
     attr(out, "pool") <- pool
     out

Modified: branches/1.15/R/specpool2vect.R
===================================================================
--- branches/1.15/R/specpool2vect.R	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/R/specpool2vect.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -1,5 +1,5 @@
 "specpool2vect" <-
-function(X, index = c("Jack.1","Jack.2", "Chao", "Boot", "Species"))
+function(X, index = c("jack1","jack2", "chao", "boot", "Species"))
 {
     pool <- attr(X, "pool")
     index <- match.arg(index)

Copied: branches/1.15/R/summary.poolaccum.R (from rev 923, pkg/vegan/R/summary.poolaccum.R)
===================================================================
--- branches/1.15/R/summary.poolaccum.R	                        (rev 0)
+++ branches/1.15/R/summary.poolaccum.R	2009-08-26 15:17:42 UTC (rev 954)
@@ -0,0 +1,24 @@
+`summary.poolaccum` <-
+    function(object, display, alpha = 0.05, ...)
+{
+    probs <- c(alpha/2, 1-alpha/2)
+    if (inherits(object, "estaccumR"))
+        dislabels <- c("S", "chao", "ace")
+    else
+        dislabels <- c("S", "chao", "jack1", "jack2", "boot")
+    disnames <- colnames(object$means[,-1])
+    names(disnames) <- dislabels
+    if (missing(display))
+        display <- dislabels
+    else
+        display <- match.arg(display, dislabels, several.ok = TRUE)
+    out <- list()
+    for (item in display) {
+        out[[item]] <- cbind(`N` = object$N,
+                             `Mean` = object$means[,disnames[item], drop=FALSE],
+                             t(apply(object[[item]], 1, quantile, probs=probs)),
+                             `Std.Dev` = apply(object[[item]], 1, sd))
+    }
+    class(out) <- "summary.poolaccum"
+    out
+}

Modified: branches/1.15/inst/ChangeLog
===================================================================
--- branches/1.15/inst/ChangeLog	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/inst/ChangeLog	2009-08-26 15:17:42 UTC (rev 954)
@@ -49,6 +49,9 @@
 
 	* merged rev685: add keyword 'datagen' in Rd where appropriate.
 
+	* added poolaccum & estaccumR with associated functions and
+	updates (rev 913 to 923).
+
 Version 1.15-3 (released June 17, 2009)
 
 	* merged revs 866 to 868: changed the way capscale displays

Modified: branches/1.15/man/specaccum.Rd
===================================================================
--- branches/1.15/man/specaccum.Rd	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/man/specaccum.Rd	2009-08-26 15:17:42 UTC (rev 954)
@@ -13,7 +13,7 @@
 }
 \usage{
 specaccum(comm, method = "exact", permutations = 100,
-          conditioned =TRUE, gamma = "Jack.1",  ...)
+          conditioned =TRUE, gamma = "jack1",  ...)
 \method{plot}{specaccum}(x, add = FALSE, ci = 2, ci.type = c("bar", "line", "polygon"), 
     col = par("fg"), ci.col = col, ci.lty = 1, xlab = "Sites", 
     ylab = x$method, ylim, ...)
@@ -136,10 +136,14 @@
   without replacement.  Further, its standard deviation does not take
   into account species correlations, and is generally too low. }
 
-\seealso{\code{\link{rarefy}} and \code{\link{renyiaccum}}.
-  Underlying graphical functions are
-  \code{\link{boxplot}}, \code{\link{matlines}}, \code{\link{segments}}
-    and \code{\link{polygon}}. }
+\seealso{\code{\link{rarefy}} and \code{\link{rrarefy}} are related
+  individual based models. Other accumulation models are
+  \code{\link{poolaccum}} for extrapoltated richenss, and
+  \code{\link{renyiaccum}} and \code{\link{tsallisaccum}} for
+  diversity indices.  Underlying graphical functions are
+  \code{\link{boxplot}}, \code{\link{matlines}},
+  \code{\link{segments}} and \code{\link{polygon}}.
+}
 \examples{
 data(BCI)
 sp1 <- specaccum(BCI)

Modified: branches/1.15/man/specpool.Rd
===================================================================
--- branches/1.15/man/specpool.Rd	2009-08-26 14:58:43 UTC (rev 953)
+++ branches/1.15/man/specpool.Rd	2009-08-26 15:17:42 UTC (rev 954)
@@ -1,10 +1,16 @@
 \name{specpool}
 \alias{specpool}
 \alias{specpool2vect}
+\alias{poolaccum}
+\alias{print.poolaccum}
+\alias{summary.poolaccum}
+\alias{plot.poolaccum}
 \alias{estimateR}
 \alias{estimateR.default}
 \alias{estimateR.matrix}
 \alias{estimateR.data.frame}
+\alias{estaccum}
+
 \title{ Extrapolated Species Richness in a Species Pool}
 \description{
   The functions estimate the extrapolated species richness in a species
@@ -15,16 +21,26 @@
 }
 \usage{
 specpool(x, pool)
-specpool2vect(X, index = c("Jack.1","Jack.2", "Chao", "Boot","Species"))
 estimateR(x, ...)
+specpool2vect(X, index = c("jack1","jack2", "chao", "boot","Species"))
+poolaccum(x, permutations = 100, minsize = 3)
+\method{summary}{poolaccum}(object, display, alpha = 0.05, ...)
+\method{plot}{poolaccum}(x, alpha = 0.05, type = c("l","g"), ...)
 }
 
 \arguments{
-  \item{x}{Data frame or matrix with species data.}
+  \item{x}{Data frame or matrix with species data or the analysis result 
+    for \code{plot} function.}
   \item{pool}{A vector giving a classification for pooling the sites in
     the species data. If missing, all sites are pooled together.}
-  \item{X}{A \code{specpool} result object.}
+  \item{X, object}{A \code{specpool} result object.}
   \item{index}{The selected index of extrapolated richness.}
+  \item{permutations}{Number of permutations of sampling order of sites.}
+  \item{minsize}{Smallest number of sampling units reported.}
+  \item{display}{Indices to be displayed.}
+  \item{alpha}{Level of quantiles shown. This proportion will be left outside
+    symmetric limits.}
+  \item{type}{Type of graph produced in \code{\link[lattice]{xyplot}}.}
   \item{...}{Other parameters (not used).}
 }
 \details{
@@ -58,6 +74,17 @@
       (1-p_i)^N}
     }
 
+  Function \code{poolaccum} is similar to \code{\link{specaccum}}, but
+  it estimates all extrapolated richness values of \code{specpool} in
+  addition to number of species for random ordering of sampling units.
+  The function has \code{summary} and \code{plot} methods. The
+  \code{summary} returns quantile envilopes of permutations
+  corresponding the given level of \code{alpha} and standard deviation
+  of permutations for each sample size. The \code{plot} function shows
+  the mean and envelope of permutations with given \code{alpha} for
+  models. The selection of models can be restricted and order changes
+  using the \code{display} argument in \code{summary} or \code{plot}.
+
     The abundance-based estimates in \code{estimateR} use counts (frequencies) of
     species in a single site. If called for a matrix or data frame, the
     function will give separate estimates for each site.  The two
@@ -105,13 +132,17 @@
     
 }
 \value{
+
   Function \code{specpool} returns a data frame with entries for
- observed richness
-  and each of the indices for each class in \code{pool} vector.  The
-  utility function \code{specpool2vect} maps the pooled values into
-  a vector giving the value of selected \code{index} for each original
-  site. Function \code{estimateR} returns the estimates and their
-  standard errors for each site. 
+  observed richness and each of the indices for each class in
+  \code{pool} vector.  The utility function \code{specpool2vect} maps
+  the pooled values into a vector giving the value of selected
+  \code{index} for each original site.  Function \code{poolaccum}
+  returns matrices of permutation results for each richness estimator,
+  the vector of sample sizes and a table of \code{means} of
+  permutations for each estimator.  Function \code{estimateR} returns
+  the estimates and their standard errors for each site.
+ 
 }
 \references{
   Chao, A. (1987). Estimating the population size for capture-recapture
@@ -137,7 +168,8 @@
   See \url{http://viceroy.eeb.uconn.edu/EstimateS} for a more complete
   (and positive) discussion and alternative software for some platforms.
 }
-\seealso{\code{\link{veiledspec}}, \code{\link{diversity}}, \code{\link{beals}}. }
+\seealso{\code{\link{veiledspec}}, \code{\link{diversity}}, \code{\link{beals}},
+ \code{\link{specaccum}}. }
 \examples{
 data(dune)
 data(dune.env)
@@ -151,6 +183,11 @@
  border="cyan3", notch=TRUE)
 par(op)
 data(BCI)
+## Accumulation model
+pool <- poolaccum(BCI)
+summary(pool, display = "chao")
+plot(pool)
+## Quantitative model
 estimateR(BCI[1:5,])
 }
 \keyword{ univar }



More information about the Vegan-commits mailing list