[Vegan-commits] r1823 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Sep 10 00:26:48 CEST 2011
Author: psolymos
Date: 2011-09-10 00:26:47 +0200 (Sat, 10 Sep 2011)
New Revision: 1823
Added:
pkg/vegan/R/clamtest.R
pkg/vegan/R/plot.clamtest.R
pkg/vegan/R/print.summary.clamtest.R
pkg/vegan/R/summary.clamtest.R
Log:
clamtest after Chazdon et al. (2011 Ecology, 92, 1332--1343)
Added: pkg/vegan/R/clamtest.R
===================================================================
--- pkg/vegan/R/clamtest.R (rev 0)
+++ pkg/vegan/R/clamtest.R 2011-09-09 22:26:47 UTC (rev 1823)
@@ -0,0 +1,112 @@
+## CLAM, reproduction of software described in Chazdon et al. 2011
+## Ecology, 92, 1332--1343
+clamtest <-
+function(comm, groups, coverage.limit = 10,
+specialization = 0.667, npoints = 20, alpha = 0.05/npoints)
+{
+ ## inital checks
+ comm <- as.matrix(comm)
+ if (NROW(comm) < 2)
+ stop("'comm' must have at least 2 rows")
+ if (nrow(comm) > 2 && missing(groups))
+ stop("'groups' is missing")
+ if (nrow(comm) == 2 && missing(groups))
+ groups <- if (is.null(rownames(comm)))
+ c("Group.1", "Group.2") else rownames(comm)
+ if (length(groups) != nrow(comm))
+ stop("length of 'groups' must equal 'nrow(comm)'")
+ if (length(unique(groups)) != 2)
+ stop("number of groups must be 2")
+ glabel <- as.character(unique(groups))
+ if (is.null(colnames(comm)))
+ colnames(comm) <- paste("Species", 1:nrow(comm), sep=".")
+ if (any(colSums(comm) <= 0))
+ stop("'comm' contains zero sum columns")
+ spp <- colnames(comm)
+ ## reproduced from Chazdon et al. 2011, Ecology 92, 1332--1343
+ S <- ncol(comm)
+ if (nrow(comm) == 2) {
+ Y <- comm[glabel[1],]
+ X <- comm[glabel[2],]
+ } else {
+ Y <- colSums(comm[which(groups==glabel[1]),])
+ X <- colSums(comm[which(groups==glabel[2]),])
+ }
+ names(X) <- names(Y) <- NULL
+ #all(ct$Total_SG == Y)
+ #all(ct$Total_OG == X)
+ m <- sum(Y)
+ n <- sum(X)
+ if (sum(Y) <= 0 || sum(X) <= 0)
+ stop("zero group totals not allowed")
+ ## check if comm contains integer, especially for singletons
+ if (any(X[X>0] < 1) || any(Y[Y>0] < 1))
+ warnings("<1 non integer values detected: analysis might not be meaningful")
+ if (abs(sum(X,Y) - sum(as.integer(X), as.integer(Y))) > 10^-6)
+ warnings("non integer values detected")
+ C1 <- 1 - sum(X==1)/n
+ C2 <- 1 - sum(Y==1)/m
+ ## this stands for other than 2/3 cases
+ uu <- specialization/(1-specialization)
+ ## critical level
+ Zp <- qnorm(alpha, lower.tail=FALSE)
+ #p_i=a
+ #pi_i=b
+ ## function to calculate test statistic from Appendix D
+ ## (Ecological Archives E092-112-A4)
+ testfun <- function(p_i, pi_i, C1, C2, n, m) {
+ C1 <- ifelse(p_i < coverage.limit, C1, 1)
+ C2 <- ifelse(pi_i < coverage.limit, C2, 1)
+ Var <- C1^2*(p_i*(1-p_i)/n) + uu^2*C2^2*(pi_i*(1-pi_i)/m)
+ C1*p_i - C2*pi_i*uu - Zp*sqrt(Var)
+ }
+ ## root finding for iso-lines (instead of itarative search)
+ rootfun <- function(pi_i, C1, C2, n, m, upper) {
+ f <- function(p_i) testfun(p_i/n, pi_i/m, C1, C2, n, m)
+ if (length(unique(sign(c(f(1), f(upper))))) > 1)
+ ceiling(uniroot(f, lower=1, upper=upper)$root) else NA
+ }
+ ## sequences for finding Xmin and Ymin values
+ Xseq <- as.integer(trunc(seq(1, max(X), len=npoints)))
+ Yseq <- as.integer(trunc(seq(1, max(Y), len=npoints)))
+ ## finding Xmin and Ymin values for Xseq and Yseq
+ Xmins <- sapply(Yseq, function(z) rootfun(z, C1, C2, n, m, upper=max(X)))
+ Ymins <- sapply(Xseq, function(z) rootfun(z, C2, C1, m, n, upper=max(Y)))
+ minval <- list(data.frame(x=Xseq[!is.na(Ymins)], y=Ymins[!is.na(Ymins)]),
+ data.frame(x=Xmins[!is.na(Xmins)], y=Yseq[!is.na(Xmins)]))
+ ## shared but too rare
+ Ymin <- Ymins[1]
+ Xmin <- Xmins[1]
+ sr <- X < Xmin & Y < Ymin
+ tmp1 <- approx(c(Xmin, 1), c(1, Ymin), xout=1:Xmin)
+ tmp2 <- approx(c(1, Ymin), c(Xmin, 1), xout=1:Ymin)
+ for (i in 1:S) {
+ if (X[i] %in% tmp1$x)
+ sr[i] <- Y[i] < tmp1$y[which(X[i]==tmp1$x)]
+ if (Y[i] %in% tmp2$x)
+ sr[i] <- X[i] < tmp2$y[which(Y[i]==tmp2$x)]
+ }
+ ## classification
+ a <- ifelse(X==0, 1, X)/n # \hat{p_i}
+ b <- ifelse(Y==0, 1, Y)/m # \hat{\pi_i}
+ specX <- !sr & testfun(a, b, C1, C2, n, m) > 0
+ specY <- !sr & testfun(b, a, C2, C1, m, n) > 0
+ gen <- !sr & !specX & !specY
+ ## crosstable
+ tmp <- ifelse(cbind(gen, specY, specX, sr), 1, 0)
+ types <- c("Generalist", paste("Specialist", glabel[1], sep="_"),
+ paste("Specialist", glabel[2], sep="_"), "Too_rare")
+ classes <- as.factor((1:4)[rowSums(tmp*col(tmp))])
+ levels(classes) <- c("Generalist", paste("Specialist", glabel[1], sep="_"),
+ paste("Specialist", glabel[2], sep="_"), "Too_rare")
+ tab <- data.frame(Species=spp, y=Y, x=X, Classes=classes)
+ colnames(tab)[2:3] <- paste("Total", glabel, sep="_")
+ rownames(tab) <- NULL
+ class(tab) <- c("clamtest","data.frame")
+ attr(tab, "settings") <- list(labels = glabel,
+ coverage.limit = coverage.limit, specialization = specialization,
+ npoints = npoints, alpha = alpha)
+ attr(tab, "minv") <- minval
+ attr(tab, "coverage") <- structure(c(C2, C1), .Names=glabel)
+ tab
+}
Added: pkg/vegan/R/plot.clamtest.R
===================================================================
--- pkg/vegan/R/plot.clamtest.R (rev 0)
+++ pkg/vegan/R/plot.clamtest.R 2011-09-09 22:26:47 UTC (rev 1823)
@@ -0,0 +1,39 @@
+plot.clamtest <- function(x, xlab, ylab, main,
+ pch=21:24, col.points=1:4, col.lines=2:4, lty=1:3,
+ position="bottomright", ...) {
+ summ <- summary(x)
+ glabel <- summ$labels
+ if (missing(main))
+ main <- "Species Classification"
+ if (missing(xlab))
+ xlab <- paste(glabel[2], "(abundance + 1)")
+ if (missing(ylab))
+ ylab <- paste(glabel[1], "(abundance + 1)")
+ Y <- x[,2]
+ X <- x[,3]
+ minval <- summ$minv
+ ## plot the dots
+ rr <- range(X+1,Y+1)
+ plot(X+1, Y+1, log = "xy", xaxt = "n", yaxt = "n",
+ col=col.points[as.integer(x$Classes)],
+ pch=pch[as.integer(x$Classes)],
+ xlab=xlab, ylab=ylab, main=main,
+ xlim=rr, ylim=rr, ...)
+ axis(1, c(1,10,100,1000,10000))
+ axis(2, c(1,10,100,1000,10000))
+ ## too rare threshold
+ Ymin <- minval[[1]][1,2]
+ Xmin <- minval[[2]][1,1]
+ lines(rep(Xmin, 2)+1, c(0, 1)+1, col=col.lines[1], lty=lty[1])
+ lines(c(0, 1)+1, rep(Ymin, 2)+1, col=col.lines[1], lty=lty[1])
+ tmp <- approx(c(Xmin, 1), c(1, Ymin))
+ lines(tmp$x+1, tmp$y+1, col=col.lines[1], lty=lty[1])
+ ## Y vs. gen threshold
+ lines(minval[[1]]+1, col=col.lines[2], lty=lty[2])
+ ## X vs. gen threshold
+ lines(minval[[2]]+1, col=col.lines[3], lty=lty[3])
+ if (!is.null(position))
+ legend(position, col=col.points, pch=pch,
+ legend=rownames(summ$summary))
+ invisible(x)
+}
Added: pkg/vegan/R/print.summary.clamtest.R
===================================================================
--- pkg/vegan/R/print.summary.clamtest.R (rev 0)
+++ pkg/vegan/R/print.summary.clamtest.R 2011-09-09 22:26:47 UTC (rev 1823)
@@ -0,0 +1,13 @@
+print.summary.clamtest <- function(x, digits=max(3, getOption("digits") - 3), ...) {
+ cat("Two Groups Species Classification Method (CLAM)\n\n")
+ cat("Specialization threshold =", x$specialization)
+ cat("\nAlpha level =", x$alpha)
+ cat("\n\nEstimated sample coverage:\n")
+ print(x$coverage, digits=digits)
+ cat("\nMinimum abundance for classification:\n")
+ print(structure(c(x$minv[[1]][1,2], x$minv[[2]][1,1]),
+ .Names=x$labels))
+ cat("\n")
+ printCoefmat(x$summary, digits=digits, ...)
+}
+
Added: pkg/vegan/R/summary.clamtest.R
===================================================================
--- pkg/vegan/R/summary.clamtest.R (rev 0)
+++ pkg/vegan/R/summary.clamtest.R 2011-09-09 22:26:47 UTC (rev 1823)
@@ -0,0 +1,7 @@
+summary.clamtest <- function(object, ...) {
+ structure(c(attr(object, "settings"),
+ list(summary=cbind(Species=table(object$Classes),
+ Proportion=table(object$Classes)/nrow(object)),
+ minv=attr(object, "minv"),
+ coverage=attr(object, "coverage"))), class="summary.clamtest")
+}
More information about the Vegan-commits
mailing list