[Vegan-commits] r2902 - in pkg/vegan: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 16 10:50:21 CEST 2014
Author: jarioksa
Date: 2014-10-16 10:50:21 +0200 (Thu, 16 Oct 2014)
New Revision: 2902
Modified:
pkg/vegan/R/clamtest.R
pkg/vegan/R/commsim.R
pkg/vegan/R/contribdiv.R
pkg/vegan/R/eventstar.R
pkg/vegan/R/indpower.R
pkg/vegan/R/make.commsim.R
pkg/vegan/R/nullmodel.R
pkg/vegan/R/permutest.betadisper.R
pkg/vegan/R/persp.tsallisaccum.R
pkg/vegan/R/plot.clamtest.R
pkg/vegan/R/plot.contribdiv.R
pkg/vegan/R/print.commsim.R
pkg/vegan/R/print.mantel.correlog.R
pkg/vegan/R/print.nullmodel.R
pkg/vegan/R/print.simmat.R
pkg/vegan/R/print.summary.clamtest.R
pkg/vegan/R/simulate.nullmodel.R
pkg/vegan/R/summary.clamtest.R
pkg/vegan/R/tsallis.R
pkg/vegan/R/tsallisaccum.R
pkg/vegan/R/update.nullmodel.R
pkg/vegan/man/clamtest.Rd
pkg/vegan/man/commsim.Rd
pkg/vegan/man/eventstar.Rd
pkg/vegan/man/mantel.correlog.Rd
pkg/vegan/man/nullmodel.Rd
pkg/vegan/man/pcnm.Rd
pkg/vegan/man/permatfull.Rd
Log:
Merge branch 'master' into r-forge-svn-local
Modified: pkg/vegan/R/clamtest.R
===================================================================
--- pkg/vegan/R/clamtest.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/clamtest.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,127 +1,127 @@
-## CLAM, reproduction of software described in Chazdon et al. 2011
-## Ecology, 92, 1332--1343
-clamtest <-
-function(comm, groups, coverage.limit = 10,
-specialization = 2/3, npoints = 20, alpha = 0.05/20)
-{
- ## 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:ncol(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))
- warning("<1 non integer values detected: analysis might not be meaningful")
- if (abs(sum(X,Y) - sum(as.integer(X), as.integer(Y))) > 10^-6)
- warning("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)
- ## coverage limit is count, not freq !!!
- testfun <- function(p_i, pi_i, C1, C2, n, m) {
- C1 <- ifelse(p_i*n < coverage.limit, C1, 1)
- C2 <- ifelse(pi_i*m < 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)))
-
- ## needed to tweak original set of rules (extreme case reported
- ## by Richard Telford failed here)
- if (all(is.na(Xmins)))
- Xmins[1] <- 1
- if (all(is.na(Ymins)))
- Ymins[1] <- 1
-
- 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
-
- ## consequence of manually setting Xmin/Ymin resolved here
- tmp1 <- if (Xmin==1)
- list(x=1, y=Xmin) else approx(c(Xmin, 1), c(1, Ymin), xout=1:Xmin)
- tmp2 <- if (Ymin==1)
- list(x=1, y=Ymin) else 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 <- factor((1:4)[rowSums(tmp*col(tmp))], levels=1:4)
- 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
-}
+## CLAM, reproduction of software described in Chazdon et al. 2011
+## Ecology, 92, 1332--1343
+clamtest <-
+function(comm, groups, coverage.limit = 10,
+specialization = 2/3, npoints = 20, alpha = 0.05/20)
+{
+ ## 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:ncol(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))
+ warning("<1 non integer values detected: analysis might not be meaningful")
+ if (abs(sum(X,Y) - sum(as.integer(X), as.integer(Y))) > 10^-6)
+ warning("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)
+ ## coverage limit is count, not freq !!!
+ testfun <- function(p_i, pi_i, C1, C2, n, m) {
+ C1 <- ifelse(p_i*n < coverage.limit, C1, 1)
+ C2 <- ifelse(pi_i*m < 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)))
+
+ ## needed to tweak original set of rules (extreme case reported
+ ## by Richard Telford failed here)
+ if (all(is.na(Xmins)))
+ Xmins[1] <- 1
+ if (all(is.na(Ymins)))
+ Ymins[1] <- 1
+
+ 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
+
+ ## consequence of manually setting Xmin/Ymin resolved here
+ tmp1 <- if (Xmin==1)
+ list(x=1, y=Xmin) else approx(c(Xmin, 1), c(1, Ymin), xout=1:Xmin)
+ tmp2 <- if (Ymin==1)
+ list(x=1, y=Ymin) else 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 <- factor((1:4)[rowSums(tmp*col(tmp))], levels=1:4)
+ 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
+}
Modified: pkg/vegan/R/commsim.R
===================================================================
--- pkg/vegan/R/commsim.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/commsim.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,24 +1,24 @@
-## this is function to create a commsim object, does some checks
-## there is a finite number of useful arguments here
-## but I added ... to allow for unforeseen algorithms,
-## or being able to reference to external objects
-commsim <-
-function(method, fun, binary, isSeq, mode)
-{
- fun <- if (!missing(fun))
- match.fun(fun) else stop("'fun' missing")
- if (any(!(names(formals(fun)) %in%
- c("x", "n", "nr", "nc", "rs", "cs", "rf", "cf", "s", "fill", "thin", "..."))))
- stop("unexpected arguments in 'fun'")
- out <- structure(list(method = if (!missing(method))
- as.character(method)[1L] else stop("'method' missing"),
- binary = if (!missing(binary))
- as.logical(binary)[1L] else stop("'binary' missing"),
- isSeq = if (!missing(isSeq))
- as.logical(isSeq)[1L] else stop("'isSeq' missing"),
- mode = if (!missing(mode))
- match.arg(as.character(mode)[1L],
- c("integer", "double")) else stop("'mode' missing"),
- fun = fun), class = "commsim")
- out
-}
+## this is function to create a commsim object, does some checks
+## there is a finite number of useful arguments here
+## but I added ... to allow for unforeseen algorithms,
+## or being able to reference to external objects
+commsim <-
+function(method, fun, binary, isSeq, mode)
+{
+ fun <- if (!missing(fun))
+ match.fun(fun) else stop("'fun' missing")
+ if (any(!(names(formals(fun)) %in%
+ c("x", "n", "nr", "nc", "rs", "cs", "rf", "cf", "s", "fill", "thin", "..."))))
+ stop("unexpected arguments in 'fun'")
+ out <- structure(list(method = if (!missing(method))
+ as.character(method)[1L] else stop("'method' missing"),
+ binary = if (!missing(binary))
+ as.logical(binary)[1L] else stop("'binary' missing"),
+ isSeq = if (!missing(isSeq))
+ as.logical(isSeq)[1L] else stop("'isSeq' missing"),
+ mode = if (!missing(mode))
+ match.arg(as.character(mode)[1L],
+ c("integer", "double")) else stop("'mode' missing"),
+ fun = fun), class = "commsim")
+ out
+}
Modified: pkg/vegan/R/contribdiv.R
===================================================================
--- pkg/vegan/R/contribdiv.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/contribdiv.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,52 +1,52 @@
-## Contribution diversity
-## Lu, H.P., H.H. Wagner and X.Y. Chen (2007).
-## A contribution diversity approach to evaluate species diversity.
-## Basic and Applied Ecology 8: 1 -12.
-`contribdiv` <-
- function(comm, index = c("richness", "simpson"), relative = FALSE,
- scaled = TRUE, drop.zero = FALSE)
-{
-
- index <- match.arg(index)
-
- x <- comm[rowSums(comm) > 0, colSums(comm) > 0]
- n <- nrow(x)
- S <- ncol(x)
-
- if (index == "richness") {
- n.i <- colSums(x > 0)
- S.k <- rowSums(x > 0)
- alpha <- S.k / n
- beta <- apply(x, 1, function(z) sum((n - n.i[z > 0]) / (n * n.i[z > 0])))
- denom <- 1
- } else {
- P.ik <- decostand(x, "total")
- P.i <- apply(P.ik, 2, function(z) sum(z) / n)
- P.i2 <- matrix(P.i, n, S, byrow=TRUE)
- alpha <- diversity(x, "simpson")
- beta <- rowSums(P.ik * (P.ik - P.i2))
- denom <- n
- }
- gamma <- alpha + beta
- D <- sum(beta) / sum(gamma)
- if (relative) {
- denom <- if (scaled)
- {denom * sum(gamma)} else 1
- alpha <- (alpha - mean(alpha)) / denom
- beta <- (beta - mean(beta)) / denom
- gamma <- (gamma - mean(gamma)) / denom
- }
- rval <- data.frame(alpha = alpha, beta = beta, gamma = gamma)
- if (!drop.zero && nrow(comm) != n) {
- nas <- rep(NA, nrow(comm))
- rval2 <- data.frame(alpha = nas, beta = nas, gamma = nas)
- rval2[rowSums(comm) > 0, ] <- rval
- rval <- rval2
- }
- attr(rval, "diff.coef") <- D
- attr(rval, "index") <- index
- attr(rval, "relative") <- relative
- attr(rval, "scaled") <- scaled
- class(rval) <- c("contribdiv", "data.frame")
- rval
-}
+## Contribution diversity
+## Lu, H.P., H.H. Wagner and X.Y. Chen (2007).
+## A contribution diversity approach to evaluate species diversity.
+## Basic and Applied Ecology 8: 1 -12.
+`contribdiv` <-
+ function(comm, index = c("richness", "simpson"), relative = FALSE,
+ scaled = TRUE, drop.zero = FALSE)
+{
+
+ index <- match.arg(index)
+
+ x <- comm[rowSums(comm) > 0, colSums(comm) > 0]
+ n <- nrow(x)
+ S <- ncol(x)
+
+ if (index == "richness") {
+ n.i <- colSums(x > 0)
+ S.k <- rowSums(x > 0)
+ alpha <- S.k / n
+ beta <- apply(x, 1, function(z) sum((n - n.i[z > 0]) / (n * n.i[z > 0])))
+ denom <- 1
+ } else {
+ P.ik <- decostand(x, "total")
+ P.i <- apply(P.ik, 2, function(z) sum(z) / n)
+ P.i2 <- matrix(P.i, n, S, byrow=TRUE)
+ alpha <- diversity(x, "simpson")
+ beta <- rowSums(P.ik * (P.ik - P.i2))
+ denom <- n
+ }
+ gamma <- alpha + beta
+ D <- sum(beta) / sum(gamma)
+ if (relative) {
+ denom <- if (scaled)
+ {denom * sum(gamma)} else 1
+ alpha <- (alpha - mean(alpha)) / denom
+ beta <- (beta - mean(beta)) / denom
+ gamma <- (gamma - mean(gamma)) / denom
+ }
+ rval <- data.frame(alpha = alpha, beta = beta, gamma = gamma)
+ if (!drop.zero && nrow(comm) != n) {
+ nas <- rep(NA, nrow(comm))
+ rval2 <- data.frame(alpha = nas, beta = nas, gamma = nas)
+ rval2[rowSums(comm) > 0, ] <- rval
+ rval <- rval2
+ }
+ attr(rval, "diff.coef") <- D
+ attr(rval, "index") <- index
+ attr(rval, "relative") <- relative
+ attr(rval, "scaled") <- scaled
+ class(rval) <- c("contribdiv", "data.frame")
+ rval
+}
Modified: pkg/vegan/R/eventstar.R
===================================================================
--- pkg/vegan/R/eventstar.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/eventstar.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,18 +1,18 @@
-eventstar <- function(x, qmax=5) {
- if (is.null(dim(x)))
- x <- matrix(x, 1, length(x))
- lossfun <- function(q, x)
- tsallis(x, scales=q, norm=TRUE)
- qstarfun <- function(x) {
- optimize(lossfun, interval=c(0, qmax), x=x)$minimum
- }
- qs <- apply(x, 1, qstarfun)
- Hs <- sapply(1:nrow(x), function(i) tsallis(x[i,],
- scales=qs[i], hill=FALSE))
- S <- rowSums(x)
- Es <- ifelse(qs==1, log(S), Hs/((S^(1-qs)-1)/(1-qs)))
- Ds <- (1-(qs-1)*Hs)^(1/(1-qs))
- out <- data.frame(qstar=qs, Estar=Es, Hstar=Hs, Dstar=Ds)
- rownames(out) <- rownames(x)
- out
-}
+eventstar <- function(x, qmax=5) {
+ if (is.null(dim(x)))
+ x <- matrix(x, 1, length(x))
+ lossfun <- function(q, x)
+ tsallis(x, scales=q, norm=TRUE)
+ qstarfun <- function(x) {
+ optimize(lossfun, interval=c(0, qmax), x=x)$minimum
+ }
+ qs <- apply(x, 1, qstarfun)
+ Hs <- sapply(1:nrow(x), function(i) tsallis(x[i,],
+ scales=qs[i], hill=FALSE))
+ S <- rowSums(x)
+ Es <- ifelse(qs==1, log(S), Hs/((S^(1-qs)-1)/(1-qs)))
+ Ds <- (1-(qs-1)*Hs)^(1/(1-qs))
+ out <- data.frame(qstar=qs, Estar=Es, Hstar=Hs, Dstar=Ds)
+ rownames(out) <- rownames(x)
+ out
+}
Modified: pkg/vegan/R/indpower.R
===================================================================
--- pkg/vegan/R/indpower.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/indpower.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,25 +1,25 @@
-indpower <-
-function(x, type=0)
-{
- x <- as.matrix(x)
- x <- ifelse(x > 0, 1, 0)
- if (NCOL(x) < 2)
- stop("provide at least 2 columns for 'x'")
- if (!(type %in% 0:2))
- stop("'type' must be in c(0, 1, 2)")
- n <- nrow(x)
- j <- crossprod(x) ## faster t(x) %*% x
- ip1 <- sweep(j, 1, diag(j), "/")
- ip2 <- 1 - sweep(-sweep(j, 2, diag(j), "-"), 1, n - diag(j), "/")
- out <- switch(as.character(type),
- "0" = sqrt(ip1 * ip2),
- "1" = ip1,
- "2" = ip2)
- cn <- if (is.null(colnames(out)))
- 1:ncol(out) else colnames(out)
- rn <- if (is.null(rownames(out)))
- 1:ncol(out) else rownames(out)
- colnames(out) <- paste("t", cn, sep=".")
- rownames(out) <- paste("i", rn, sep=".")
- out
-}
+indpower <-
+function(x, type=0)
+{
+ x <- as.matrix(x)
+ x <- ifelse(x > 0, 1, 0)
+ if (NCOL(x) < 2)
+ stop("provide at least 2 columns for 'x'")
+ if (!(type %in% 0:2))
+ stop("'type' must be in c(0, 1, 2)")
+ n <- nrow(x)
+ j <- crossprod(x) ## faster t(x) %*% x
+ ip1 <- sweep(j, 1, diag(j), "/")
+ ip2 <- 1 - sweep(-sweep(j, 2, diag(j), "-"), 1, n - diag(j), "/")
+ out <- switch(as.character(type),
+ "0" = sqrt(ip1 * ip2),
+ "1" = ip1,
+ "2" = ip2)
+ cn <- if (is.null(colnames(out)))
+ 1:ncol(out) else colnames(out)
+ rn <- if (is.null(rownames(out)))
+ 1:ncol(out) else rownames(out)
+ colnames(out) <- paste("t", cn, sep=".")
+ rownames(out) <- paste("i", rn, sep=".")
+ out
+}
Modified: pkg/vegan/R/make.commsim.R
===================================================================
--- pkg/vegan/R/make.commsim.R 2014-10-16 08:50:18 UTC (rev 2901)
+++ pkg/vegan/R/make.commsim.R 2014-10-16 08:50:21 UTC (rev 2902)
@@ -1,436 +1,436 @@
-## this lists all known algos in vegan and more
-## if method is commsim object, it is returned
-## if it is character, switch returns the right one, else stop with error
-## so it can be used instead of match.arg(method) in other functions
-## NOTE: very very long -- but it can be a central repository of algos
-## NOTE 2: storage mode coercions are avoided here
-## (with no apparent effect on speed), it should be
-## handled by nullmodel and commsim characteristics
-make.commsim <-
-function(method)
-{
- algos <- list(
- "r00" = commsim(method="r00", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- matrix(0L, nr * nc, n)
- for (k in seq_len(n))
- out[sample.int(nr * nc, s), k] <- 1L
- dim(out) <- c(nr, nc, n)
- out
- }),
- "c0" = commsim(method="c0", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- J <- seq_len(nc)
- for (k in seq_len(n))
- for (j in J)
- out[sample.int(nr, cs[j]), j, k] <- 1L
- out
- }),
- "r0" = commsim(method="r0", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- I <- seq_len(nr)
- for (k in seq_len(n))
- for (i in I)
- out[i, sample.int(nc, rs[i]), k] <- 1L
- out
- }),
- "r0_old" = commsim(method="r0_old", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- I <- seq_len(nr)
- p <- rep(1, nc)
- for (k in seq_len(n))
- for (i in I)
- out[i, sample.int(nc, rs[i], prob = p), k] <- 1L
- out
- }),
- "r1" = commsim(method="r1", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- I <- seq_len(nr)
- storage.mode(cs) <- "double"
- for (k in seq_len(n))
- for (i in I)
- out[i, sample.int(nc, rs[i], prob=cs), k] <- 1L
- out
- }),
- "r2" = commsim(method="r2", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- p <- cs * cs
- I <- seq_len(nr)
- for (k in seq_len(n))
- for (i in I)
- out[i, sample.int(nc, rs[i], prob=p), k] <- 1L
- out
- }),
- "quasiswap" = commsim(method="quasiswap", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(unlist(r2dtable(n, rs, cs)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n))
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- out
- }),
- "swap" = commsim(method="swap", binary=TRUE, isSeq=TRUE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- out[,,1] <- .C("swap",
- m = x, nr, nc, thin, PACKAGE = "vegan")$m
- for (k in seq_len(n-1))
- out[,,k+1] <- .C("swap",
- m = out[,,k], nr, nc, thin,
- PACKAGE = "vegan")$m
- out
- }),
- "tswap" = commsim(method="tswap", binary=TRUE, isSeq=TRUE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- out[,,1] <- .C("trialswap",
- m = x, nr, nc, thin, PACKAGE = "vegan")$m
- for (k in seq_len(n-1))
- out[,,k+1] <- .C("trialswap",
- m = out[,,k], nr, nc, thin, PACKAGE = "vegan")$m
- out
- }),
- "backtrack" = commsim(method="backtrack", binary=TRUE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, rs, cs, rf, cf, s, fill, thin) {
- btrfun <- function() {
- all <- matrix(as.integer(1:(nr * nc)), nrow = nr, ncol = nc)
- out <- matrix(0L, nrow = nr, ncol = nc)
- free <- matrix(as.integer(1:(nr * nc)), nrow = nr)
- icount <- integer(length(rs))
- jcount <- integer(length(cs))
- prob <- outer(rs, cs, "*")
- ij <- sample(free, prob = prob)
- i <- (ij - 1)%%nr + 1
- j <- (ij - 1)%/%nr + 1
- for (k in 1:length(ij)) {
- if (icount[i[k]] < rs[i[k]] && jcount[j[k]] < cs[j[k]]) {
- out[ij[k]] <- 1L
- icount[i[k]] <- icount[i[k]] + 1L
- jcount[j[k]] <- jcount[j[k]] + 1L
- }
- }
- ndrop <- 1
- for (i in 1:10000) {
- oldout <- out
- oldn <- sum(out)
- drop <- sample(all[out == 1L], ndrop)
- out[drop] <- 0L
- candi <- outer(rowSums(out) < rs, colSums(out) < cs, "&") & out == 0L
- while (sum(candi) > 0) {
- if (sum(candi) > 1)
- ij <- sample(all[candi], 1)
- else ij <- all[candi]
- out[ij] <- 1L
- candi <- outer(rowSums(out) < rs, colSums(out) < cs, "&") & out == 0
- }
- if (sum(out) >= fill)
- break
- if (oldn >= sum(out))
- ndrop <- min(ndrop + 1, 4)
- else ndrop <- 1
- if (oldn > sum(out))
- out <- oldout
- }
- out
- }
- out <- array(0L, c(nr, nc, n))
- for (k in seq_len(n))
- out[, , k] <- btrfun()
- out
- }),
- "r2dtable" = commsim(method="r2dtable", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(unlist(r2dtable(n, rs, cs)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- out
- }),
- "swap_count" = commsim(method="swap_count", binary=FALSE, isSeq=TRUE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(0L, c(nr, nc, n))
- out[,,1] <- .C("swapcount",
- m = x, nr, nc, thin, PACKAGE = "vegan")$m
- for (k in seq_len(n-1))
- out[,,k+1] <- .C("swapcount",
- m = out[,,k], nr, nc, thin, PACKAGE = "vegan")$m
- out
- }),
- "quasiswap_count" = commsim(method="quasiswap_count", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(unlist(r2dtable(n, rs, cs)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n))
- out[,,k] <- .C("rswapcount",
- m = out[,,k], nr, nc, fill, PACKAGE = "vegan")$m
- out
- }),
- "swsh_samp" = commsim(method="swsh_samp", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- nz <- as.integer(x[x > 0])
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- out[,,k][out[,,k] > 0] <- sample(nz) # we assume that length(nz)>1
- }
- out
- }),
- "swsh_both" = commsim(method="swsh_both", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- indshuffle <- function(x) {
- drop(rmultinom(1, sum(x), rep(1, length(x))))
- }
- nz <- as.integer(x[x > 0])
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- out[,,k][out[,,k] > 0] <- indshuffle(nz - 1L) + 1L # we assume that length(nz)>1
- }
- out
- }),
- "swsh_samp_r" = commsim(method="swsh_samp_r", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- I <- seq_len(nr)
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- for (i in I) {
- nz <- as.integer(x[i,][x[i,] > 0])
- if (length(nz) == 1)
- out[i,,k][out[i,,k] > 0] <- nz
- if (length(nz) > 1)
- out[i,,k][out[i,,k] > 0] <- sample(nz)
- }
- }
- out
- }),
- "swsh_samp_c" = commsim(method="swsh_samp_c", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- J <- seq_len(nc)
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- for (j in J) {
- nz <- as.integer(x[,j][x[,j] > 0])
- if (length(nz) == 1)
- out[,j,k][out[,j,k] > 0] <- nz
- if (length(nz) > 1)
- out[,j,k][out[,j,k] > 0] <- sample(nz)
- }
- }
- out
- }),
- "swsh_both_r" = commsim(method="swsh_both_r", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- indshuffle <- function(x) {
- drop(rmultinom(1, sum(x), rep(1, length(x))))
- }
- I <- seq_len(nr)
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- for (i in I) {
- nz <- as.integer(x[i,][x[i,] > 0])
- if (length(nz) == 1)
- out[i,,k][out[i,,k] > 0] <- nz
- if (length(nz) > 1)
- out[i,,k][out[i,,k] > 0] <- indshuffle(nz - 1L) + 1L
- }
- }
- out
- }),
- "swsh_both_c" = commsim(method="swsh_both_c", binary=FALSE, isSeq=FALSE,
- mode="integer",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- indshuffle <- function(x) {
- drop(rmultinom(1, sum(x), rep(1, length(x))))
- }
- J <- seq_len(nc)
- out <- array(unlist(r2dtable(fill, rf, cf)), c(nr, nc, n))
- storage.mode(out) <- "integer"
- for (k in seq_len(n)) {
- out[,,k] <- .C("quasiswap",
- m = out[,,k], nr, nc, PACKAGE = "vegan")$m
- for (j in J) {
- nz <- as.integer(x[,j][x[,j] > 0])
- if (length(nz) == 1)
- out[,j,k][out[,j,k] > 0] <- nz
- if (length(nz) > 1)
- out[,j,k][out[,j,k] > 0] <- indshuffle(nz - 1L) + 1L
- }
- }
- out
- }),
- "abuswap_r" = commsim(method="abuswap_r", binary=FALSE, isSeq=TRUE,
- mode="double",
- fun=function(x, n, nr, nc, cs, rs, rf, cf, s, fill, thin) {
- out <- array(0, c(nr, nc, n))
- out[,,1] <- .C("abuswap",
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vegan -r 2902
More information about the Vegan-commits
mailing list