[Vegan-commits] r2151 - in branches/2.0: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 6 15:55:16 CEST 2012
Author: jarioksa
Date: 2012-05-06 15:55:15 +0200 (Sun, 06 May 2012)
New Revision: 2151
Added:
branches/2.0/R/adipart.default.R
branches/2.0/R/adipart.formula.R
branches/2.0/R/hiersimu.default.R
branches/2.0/R/hiersimu.formula.R
branches/2.0/R/multipart.default.R
branches/2.0/R/multipart.formula.R
Modified:
branches/2.0/NAMESPACE
branches/2.0/R/adipart.R
branches/2.0/R/hiersimu.R
branches/2.0/R/multipart.R
branches/2.0/inst/ChangeLog
branches/2.0/inst/NEWS.Rd
branches/2.0/man/adipart.Rd
branches/2.0/man/multipart.Rd
Log:
merge r2143: formula and default method for adipart, hiersimu and multipart and easier trivial analysis of alpha, beta and gamma diversities
Modified: branches/2.0/NAMESPACE
===================================================================
--- branches/2.0/NAMESPACE 2012-05-02 05:55:40 UTC (rev 2150)
+++ branches/2.0/NAMESPACE 2012-05-06 13:55:15 UTC (rev 2151)
@@ -68,6 +68,9 @@
if (getRversion() < "2.13.0") {
importFrom(permute, nobs)
}
+# adipart: vegan
+S3method(adipart, default)
+S3method(adipart, formula)
# AIC: stats
S3method(AIC, radfit)
# RsquareAdj: vegan
@@ -168,6 +171,9 @@
S3method(goodness, rda)
# head: utils
S3method(head, summary.cca)
+# hiersimu: vegan
+S3method(hiersimu, default)
+S3method(hiersimu, formula)
# identify: graphics
S3method(identify, ordiplot)
# lines: graphics
@@ -181,6 +187,9 @@
# model.frame, model.matrix: stats
S3method(model.frame, cca)
S3method(model.matrix, cca)
+# multipart: vegan
+S3method(multipart, default)
+S3method(multipart, formula)
# nobs: stats
S3method(nobs, CCorA)
S3method(nobs, adonis)
Modified: branches/2.0/R/adipart.R
===================================================================
--- branches/2.0/R/adipart.R 2012-05-02 05:55:40 UTC (rev 2150)
+++ branches/2.0/R/adipart.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -1,114 +1,6 @@
adipart <-
-function(formula, data, index=c("richness", "shannon", "simpson"),
- weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
+function (...)
{
- ## evaluate formula
- lhs <- formula[[2]]
- if (missing(data))
- data <- parent.frame()
- lhs <- as.matrix(eval(lhs, data))
- formula[[2]] <- NULL
- rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
- tlab <- attr(attr(rhs, "terms"), "term.labels")
- 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]]
- int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
- if (!identical(noint, int))
- stop("interactions are not allowed in formula")
- if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
- stop("all right hand side variables in formula must be factors")
- l1 <- sapply(rhs, function(z) length(unique(z)))
- if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
- stop("number of levels are inapropriate, check sequence")
- rval <- list()
- rval[[1]] <- as.factor(rhs[,nlevs])
- rval[[1]] <- rval[[1]][drop = TRUE]
- nCol <- nlevs - 1
- for (i in 2:nlevs) {
- rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
- nCol <- nCol - 1
- }
- rval <- as.data.frame(rval[rev(1:length(rval))])
- l2 <- sapply(rval, function(z) length(unique(z)))
- if (any(l1 != l2))
- warning("levels are not perfectly nested")
-
- ## aggregate response matrix
- fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
- TRUE else FALSE
- ftmp <- vector("list", nlevs)
- for (i in 1:nlevs) {
- ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
- }
-
- ## is there a method/burnin/thin in ... ?
- method <- if (is.null(list(...)$method))
- "r2dtable" else list(...)$method
- burnin <- if (is.null(list(...)$burnin))
- 0 else list(...)$burnin
- thin <- if (is.null(list(...)$thin))
- 1 else list(...)$thin
- base <- if (is.null(list(...)$base))
- exp(1) else list(...)$base
-
- ## evaluate other arguments
- index <- match.arg(index)
- weights <- match.arg(weights)
- switch(index,
- "richness" = {
- divfun <- function(x) apply(x > 0, 1, sum)},
- "shannon" = {
- divfun <- function(x) diversity(x, index = "shannon", MARGIN = 1, base=base)},
- "simpson" = {
- divfun <- function(x) diversity(x, index = "simpson", MARGIN = 1)})
-
- ## this is the function passed to oecosimu
- wdivfun <- function(x) {
- ## matrix sum *can* change in oecosimu (but default is constant sumMatr)
- sumMatr <- sum(x)
- if (fullgamma) {
- tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
- } else {
- tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- }
- ## weights will change in oecosimu thus need to be recalculated
- if (weights == "prop")
- wt <- lapply(1:nlevs, function(i) apply(tmp[[i]], 1, function(z) sum(z) / sumMatr))
- else wt <- lapply(1:nlevs, function(i) rep(1 / NROW(tmp[[i]]), NROW(tmp[[i]])))
- a <- sapply(1:nlevs, function(i) sum(divfun(tmp[[i]]) * wt[[i]]))
- if (relative)
- a <- a / a[length(a)]
- b <- sapply(2:nlevs, function(i) a[i] - a[(i-1)])
- c(a, b)
- }
- if (nsimul > 0) {
- sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
- burnin=burnin, thin=thin)
- } 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
- attr(sim, "call") <- match.call()
- attr(sim, "index") <- index
- attr(sim, "weights") <- weights
- attr(sim, "n.levels") <- nlevs
- attr(sim, "terms") <- tlab
- attr(sim, "model") <- rhs
- class(sim) <- c("adipart", "list")
- sim
+ UseMethod("adipart")
}
Copied: branches/2.0/R/adipart.default.R (from rev 2143, pkg/vegan/R/adipart.default.R)
===================================================================
--- branches/2.0/R/adipart.default.R (rev 0)
+++ branches/2.0/R/adipart.default.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -0,0 +1,109 @@
+adipart.default <-
+function(y, x, index=c("richness", "shannon", "simpson"),
+ weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
+{
+ ## evaluate formula
+ lhs <- as.matrix(y)
+ if (missing(x))
+ x <- cbind(level_1=seq_len(nrow(lhs)),
+ leve_2=rep(1, nrow(lhs)))
+ rhs <- data.frame(x)
+ rhs[] <- lapply(rhs, as.factor)
+ rhs[] <- lapply(rhs, droplevels)
+ nlevs <- ncol(rhs)
+ 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")
+ if (is.null(colnames(rhs)))
+ colnames(rhs) <- paste("level", 1:nlevs, sep="_")
+ tlab <- colnames(rhs)
+
+ ## part check proper design of the model frame
+ l1 <- sapply(rhs, function(z) length(unique(z)))
+ if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
+ stop("number of levels are inapropriate, check sequence")
+ rval <- list()
+ rval[[1]] <- rhs[,nlevs]
+ nCol <- nlevs - 1
+ for (i in 2:nlevs) {
+ rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
+ nCol <- nCol - 1
+ }
+ rval <- as.data.frame(rval[rev(1:length(rval))])
+ l2 <- sapply(rval, function(z) length(unique(z)))
+ if (any(l1 != l2))
+ warning("levels are not perfectly nested")
+
+ ## aggregate response matrix
+ fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
+ TRUE else FALSE
+ ftmp <- vector("list", nlevs)
+ for (i in 1:nlevs) {
+ ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
+ }
+
+ ## is there a method/burnin/thin in ... ?
+ method <- if (is.null(list(...)$method))
+ "r2dtable" else list(...)$method
+ burnin <- if (is.null(list(...)$burnin))
+ 0 else list(...)$burnin
+ thin <- if (is.null(list(...)$thin))
+ 1 else list(...)$thin
+ base <- if (is.null(list(...)$base))
+ exp(1) else list(...)$base
+
+ ## evaluate other arguments
+ index <- match.arg(index)
+ weights <- match.arg(weights)
+ switch(index,
+ "richness" = {
+ divfun <- function(x) apply(x > 0, 1, sum)},
+ "shannon" = {
+ divfun <- function(x) diversity(x, index = "shannon", MARGIN = 1, base=base)},
+ "simpson" = {
+ divfun <- function(x) diversity(x, index = "simpson", MARGIN = 1)})
+
+ ## this is the function passed to oecosimu
+ wdivfun <- function(x) {
+ ## matrix sum *can* change in oecosimu (but default is constant sumMatr)
+ sumMatr <- sum(x)
+ if (fullgamma) {
+ tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
+ } else {
+ tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ }
+ ## weights will change in oecosimu thus need to be recalculated
+ if (weights == "prop")
+ wt <- lapply(1:nlevs, function(i) apply(tmp[[i]], 1, function(z) sum(z) / sumMatr))
+ else wt <- lapply(1:nlevs, function(i) rep(1 / NROW(tmp[[i]]), NROW(tmp[[i]])))
+ a <- sapply(1:nlevs, function(i) sum(divfun(tmp[[i]]) * wt[[i]]))
+ if (relative)
+ a <- a / a[length(a)]
+ b <- sapply(2:nlevs, function(i) a[i] - a[(i-1)])
+ c(a, b)
+ }
+ if (nsimul > 0) {
+ sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
+ burnin=burnin, thin=thin)
+ } 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
+ attr(sim, "call") <- match.call()
+ attr(sim, "index") <- index
+ attr(sim, "weights") <- weights
+ attr(sim, "n.levels") <- nlevs
+ attr(sim, "terms") <- tlab
+ attr(sim, "model") <- rhs
+ class(sim) <- c("adipart", "list")
+ sim
+}
Copied: branches/2.0/R/adipart.formula.R (from rev 2143, pkg/vegan/R/adipart.formula.R)
===================================================================
--- branches/2.0/R/adipart.formula.R (rev 0)
+++ branches/2.0/R/adipart.formula.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -0,0 +1,113 @@
+adipart.formula <-
+function(formula, data, index=c("richness", "shannon", "simpson"),
+ weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
+{
+ ## evaluate formula
+ lhs <- formula[[2]]
+ if (missing(data))
+ data <- parent.frame()
+ lhs <- as.matrix(eval(lhs, data))
+ formula[[2]] <- NULL
+ rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
+ tlab <- attr(attr(rhs, "terms"), "term.labels")
+ 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]]
+ int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
+ if (!identical(noint, int))
+ stop("interactions are not allowed in formula")
+ if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
+ stop("all right hand side variables in formula must be factors")
+ l1 <- sapply(rhs, function(z) length(unique(z)))
+ if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
+ stop("number of levels are inapropriate, check sequence")
+ rval <- list()
+ rval[[1]] <- as.factor(rhs[,nlevs])
+ rval[[1]] <- rval[[1]][drop = TRUE]
+ nCol <- nlevs - 1
+ for (i in 2:nlevs) {
+ rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
+ nCol <- nCol - 1
+ }
+ rval <- as.data.frame(rval[rev(1:length(rval))])
+ l2 <- sapply(rval, function(z) length(unique(z)))
+ if (any(l1 != l2))
+ warning("levels are not perfectly nested")
+
+ ## aggregate response matrix
+ fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
+ TRUE else FALSE
+ ftmp <- vector("list", nlevs)
+ for (i in 1:nlevs) {
+ ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
+ }
+
+ ## is there a method/burnin/thin in ... ?
+ method <- if (is.null(list(...)$method))
+ "r2dtable" else list(...)$method
+ burnin <- if (is.null(list(...)$burnin))
+ 0 else list(...)$burnin
+ thin <- if (is.null(list(...)$thin))
+ 1 else list(...)$thin
+ base <- if (is.null(list(...)$base))
+ exp(1) else list(...)$base
+
+ ## evaluate other arguments
+ index <- match.arg(index)
+ weights <- match.arg(weights)
+ switch(index,
+ "richness" = {
+ divfun <- function(x) apply(x > 0, 1, sum)},
+ "shannon" = {
+ divfun <- function(x) diversity(x, index = "shannon", MARGIN = 1, base=base)},
+ "simpson" = {
+ divfun <- function(x) diversity(x, index = "simpson", MARGIN = 1)})
+
+ ## this is the function passed to oecosimu
+ wdivfun <- function(x) {
+ ## matrix sum *can* change in oecosimu (but default is constant sumMatr)
+ sumMatr <- sum(x)
+ if (fullgamma) {
+ tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
+ } else {
+ tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ }
+ ## weights will change in oecosimu thus need to be recalculated
+ if (weights == "prop")
+ wt <- lapply(1:nlevs, function(i) apply(tmp[[i]], 1, function(z) sum(z) / sumMatr))
+ else wt <- lapply(1:nlevs, function(i) rep(1 / NROW(tmp[[i]]), NROW(tmp[[i]])))
+ a <- sapply(1:nlevs, function(i) sum(divfun(tmp[[i]]) * wt[[i]]))
+ if (relative)
+ a <- a / a[length(a)]
+ b <- sapply(2:nlevs, function(i) a[i] - a[(i-1)])
+ c(a, b)
+ }
+ if (nsimul > 0) {
+ sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
+ burnin=burnin, thin=thin)
+ } 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
+ attr(sim, "call") <- match.call()
+ attr(sim, "index") <- index
+ attr(sim, "weights") <- weights
+ attr(sim, "n.levels") <- nlevs
+ attr(sim, "terms") <- tlab
+ attr(sim, "model") <- rhs
+ class(sim) <- c("adipart", "list")
+ sim
+}
Modified: branches/2.0/R/hiersimu.R
===================================================================
--- branches/2.0/R/hiersimu.R 2012-05-02 05:55:40 UTC (rev 2150)
+++ branches/2.0/R/hiersimu.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -1,96 +1,5 @@
hiersimu <-
-function(formula, data, FUN, location = c("mean", "median"),
-relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
+function (...)
{
- ## evaluate formula
- lhs <- formula[[2]]
- if (missing(data))
- data <- parent.frame()
- lhs <- as.matrix(eval(lhs, data))
- formula[[2]] <- NULL
- rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
- tlab <- attr(attr(rhs, "terms"), "term.labels")
- nlevs <- length(tlab)
-
- ## part check proper design of the model frame
- noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
- int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
- if (!identical(noint, int))
- stop("interactions are not allowed in formula")
- if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
- stop("all right hand side variables in formula must be factors")
- l1 <- sapply(rhs, function(z) length(unique(z)))
- if (nlevs > 1 && !any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
- stop("number of levels are inapropriate, check sequence")
- rval <- list()
- rval[[1]] <- as.factor(rhs[,nlevs])
- rval[[1]] <- rval[[1]][drop = TRUE]
- if (nlevs > 1) {
- nCol <- nlevs - 1
- for (i in 2:nlevs) {
- rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
- nCol <- nCol - 1
- }
- }
- rval <- as.data.frame(rval[rev(1:length(rval))])
- l2 <- sapply(rval, function(z) length(unique(z)))
- if (any(l1 != l2))
- warning("levels are not perfectly nested")
-
- ## aggregate response matrix
- fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
- TRUE else FALSE
- if (fullgamma && drop.highest)
- nlevs <- nlevs - 1
- if (nlevs == 1 && relative)
- stop("'relative=FALSE' makes no sense with 1 level")
- ftmp <- vector("list", nlevs)
- for (i in 1:nlevs) {
- ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
- }
-
- ## is there a method/burnin/thin in ... ?
- method <- if (is.null(list(...)$method))
- "r2dtable" else list(...)$method
- burnin <- if (is.null(list(...)$burnin))
- 0 else list(...)$burnin
- thin <- if (is.null(list(...)$thin))
- 1 else list(...)$thin
-
- ## evaluate other arguments
- if (!is.function(FUN))
- stop("'FUN' must be a function")
- location <- match.arg(location)
- aggrFUN <- switch(location,
- "mean" = mean,
- "median" = median)
-
- ## this is the function passed to oecosimu
- evalFUN <- function(x) {
- if (fullgamma && !drop.highest) {
- tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- tmp[[nlevs]] <- matrix(colSums(lhs), nrow = 1, ncol = ncol(lhs))
- } else {
- tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- }
- a <- sapply(1:nlevs, function(i) aggrFUN(FUN(tmp[[i]]))) # dots removed from FUN
- if (relative)
- a <- a / a[length(a)]
- a
- }
-
- ## processing oecosimu results
- sim <- oecosimu(lhs, evalFUN, method = method, nsimul=nsimul,
- burnin=burnin, thin=thin)
-# nam <- paste("level", 1:nlevs, sep=".")
-# names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
- names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- tlab[1:nlevs]
- attr(sim, "call") <- match.call()
- attr(sim, "FUN") <- FUN
- attr(sim, "location") <- location
- attr(sim, "n.levels") <- nlevs
- attr(sim, "terms") <- tlab
- attr(sim, "model") <- rhs
- class(sim) <- c("hiersimu", "list")
- sim
+ UseMethod("hiersimu")
}
Copied: branches/2.0/R/hiersimu.default.R (from rev 2143, pkg/vegan/R/hiersimu.default.R)
===================================================================
--- branches/2.0/R/hiersimu.default.R (rev 0)
+++ branches/2.0/R/hiersimu.default.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -0,0 +1,93 @@
+hiersimu.default <-
+function(y, x, FUN, location = c("mean", "median"),
+relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
+{
+ ## evaluate formula
+ lhs <- as.matrix(y)
+ if (missing(x))
+ x <- cbind(level_1=seq_len(nrow(lhs)),
+ leve_2=rep(1, nrow(lhs)))
+ rhs <- data.frame(x)
+ rhs[] <- lapply(rhs, as.factor)
+ rhs[] <- lapply(rhs, droplevels)
+ nlevs <- ncol(rhs)
+ if (is.null(colnames(rhs)))
+ colnames(rhs) <- paste("level", 1:nlevs, sep="_")
+ tlab <- colnames(rhs)
+
+ ## part check proper design of the model frame
+ l1 <- sapply(rhs, function(z) length(unique(z)))
+ if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
+ stop("number of levels are inapropriate, check sequence")
+ rval <- list()
+ rval[[1]] <- rhs[,nlevs]
+ nCol <- nlevs - 1
+ if (nlevs > 1) {
+ nCol <- nlevs - 1
+ for (i in 2:nlevs) {
+ rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
+ nCol <- nCol - 1
+ }
+ }
+ rval <- as.data.frame(rval[rev(1:length(rval))])
+ l2 <- sapply(rval, function(z) length(unique(z)))
+ if (any(l1 != l2))
+ warning("levels are not perfectly nested")
+
+ ## aggregate response matrix
+ fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
+ TRUE else FALSE
+ if (fullgamma && drop.highest)
+ nlevs <- nlevs - 1
+ if (nlevs == 1 && relative)
+ stop("'relative=FALSE' makes no sense with 1 level")
+ ftmp <- vector("list", nlevs)
+ for (i in 1:nlevs) {
+ ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
+ }
+
+ ## is there a method/burnin/thin in ... ?
+ method <- if (is.null(list(...)$method))
+ "r2dtable" else list(...)$method
+ burnin <- if (is.null(list(...)$burnin))
+ 0 else list(...)$burnin
+ thin <- if (is.null(list(...)$thin))
+ 1 else list(...)$thin
+
+ ## evaluate other arguments
+ if (!is.function(FUN))
+ stop("'FUN' must be a function")
+ location <- match.arg(location)
+ aggrFUN <- switch(location,
+ "mean" = mean,
+ "median" = median)
+
+ ## this is the function passed to oecosimu
+ evalFUN <- function(x) {
+ if (fullgamma && !drop.highest) {
+ tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ tmp[[nlevs]] <- matrix(colSums(lhs), nrow = 1, ncol = ncol(lhs))
+ } else {
+ tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ }
+ a <- sapply(1:nlevs, function(i) aggrFUN(FUN(tmp[[i]]))) # dots removed from FUN
+ if (relative)
+ a <- a / a[length(a)]
+ a
+ }
+
+ ## processing oecosimu results
+ sim <- oecosimu(lhs, evalFUN, method = method, nsimul=nsimul,
+ burnin=burnin, thin=thin)
+# nam <- paste("level", 1:nlevs, sep=".")
+# names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
+ names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- tlab[1:nlevs]
+ attr(sim, "call") <- match.call()
+ attr(sim, "FUN") <- FUN
+ attr(sim, "location") <- location
+ attr(sim, "n.levels") <- nlevs
+ attr(sim, "terms") <- tlab
+ attr(sim, "model") <- rhs
+ class(sim) <- c("hiersimu", "list")
+ sim
+}
Copied: branches/2.0/R/hiersimu.formula.R (from rev 2143, pkg/vegan/R/hiersimu.formula.R)
===================================================================
--- branches/2.0/R/hiersimu.formula.R (rev 0)
+++ branches/2.0/R/hiersimu.formula.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -0,0 +1,96 @@
+hiersimu.formula <-
+function(formula, data, FUN, location = c("mean", "median"),
+relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
+{
+ ## evaluate formula
+ lhs <- formula[[2]]
+ if (missing(data))
+ data <- parent.frame()
+ lhs <- as.matrix(eval(lhs, data))
+ formula[[2]] <- NULL
+ rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
+ tlab <- attr(attr(rhs, "terms"), "term.labels")
+ nlevs <- length(tlab)
+
+ ## part check proper design of the model frame
+ noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
+ int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
+ if (!identical(noint, int))
+ stop("interactions are not allowed in formula")
+ if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
+ stop("all right hand side variables in formula must be factors")
+ l1 <- sapply(rhs, function(z) length(unique(z)))
+ if (nlevs > 1 && !any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
+ stop("number of levels are inapropriate, check sequence")
+ rval <- list()
+ rval[[1]] <- as.factor(rhs[,nlevs])
+ rval[[1]] <- rval[[1]][drop = TRUE]
+ if (nlevs > 1) {
+ nCol <- nlevs - 1
+ for (i in 2:nlevs) {
+ rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
+ nCol <- nCol - 1
+ }
+ }
+ rval <- as.data.frame(rval[rev(1:length(rval))])
+ l2 <- sapply(rval, function(z) length(unique(z)))
+ if (any(l1 != l2))
+ warning("levels are not perfectly nested")
+
+ ## aggregate response matrix
+ fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
+ TRUE else FALSE
+ if (fullgamma && drop.highest)
+ nlevs <- nlevs - 1
+ if (nlevs == 1 && relative)
+ stop("'relative=FALSE' makes no sense with 1 level")
+ ftmp <- vector("list", nlevs)
+ for (i in 1:nlevs) {
+ ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
+ }
+
+ ## is there a method/burnin/thin in ... ?
+ method <- if (is.null(list(...)$method))
+ "r2dtable" else list(...)$method
+ burnin <- if (is.null(list(...)$burnin))
+ 0 else list(...)$burnin
+ thin <- if (is.null(list(...)$thin))
+ 1 else list(...)$thin
+
+ ## evaluate other arguments
+ if (!is.function(FUN))
+ stop("'FUN' must be a function")
+ location <- match.arg(location)
+ aggrFUN <- switch(location,
+ "mean" = mean,
+ "median" = median)
+
+ ## this is the function passed to oecosimu
+ evalFUN <- function(x) {
+ if (fullgamma && !drop.highest) {
+ tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ tmp[[nlevs]] <- matrix(colSums(lhs), nrow = 1, ncol = ncol(lhs))
+ } else {
+ tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
+ }
+ a <- sapply(1:nlevs, function(i) aggrFUN(FUN(tmp[[i]]))) # dots removed from FUN
+ if (relative)
+ a <- a / a[length(a)]
+ a
+ }
+
+ ## processing oecosimu results
+ sim <- oecosimu(lhs, evalFUN, method = method, nsimul=nsimul,
+ burnin=burnin, thin=thin)
+# nam <- paste("level", 1:nlevs, sep=".")
+# names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
+ names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- tlab[1:nlevs]
+ attr(sim, "call") <- match.call()
+ attr(sim, "FUN") <- FUN
+ attr(sim, "location") <- location
+ attr(sim, "n.levels") <- nlevs
+ attr(sim, "terms") <- tlab
+ attr(sim, "model") <- rhs
+ class(sim) <- c("hiersimu", "list")
+ sim
+}
Modified: branches/2.0/R/multipart.R
===================================================================
--- branches/2.0/R/multipart.R 2012-05-02 05:55:40 UTC (rev 2150)
+++ branches/2.0/R/multipart.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -1,139 +1,5 @@
multipart <-
-function(formula, data, index=c("renyi", "tsallis"), scales = 1,
- global = FALSE, relative = FALSE, nsimul=99, ...)
+function (...)
{
- if (length(scales) > 1)
- stop("length of 'scales' must be 1")
- ## evaluate formula
- lhs <- formula[[2]]
- if (missing(data))
- data <- parent.frame()
- lhs <- as.matrix(eval(lhs, data))
- formula[[2]] <- NULL
- rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
- tlab <- attr(attr(rhs, "terms"), "term.labels")
- 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]]
- int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
- if (!identical(noint, int))
- stop("interactions are not allowed in formula")
- if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
- stop("all right hand side variables in formula must be factors")
- l1 <- sapply(rhs, function(z) length(unique(z)))
- if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
- stop("number of levels are inapropriate, check sequence")
- rval <- list()
- rval[[1]] <- as.factor(rhs[,nlevs])
- rval[[1]] <- rval[[1]][drop = TRUE]
- nCol <- nlevs - 1
- for (i in 2:nlevs) {
- rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
- nCol <- nCol - 1
- }
- rval <- as.data.frame(rval[rev(1:length(rval))])
- l2 <- sapply(rval, function(z) length(unique(z)))
- if (any(l1 != l2))
- warning("levels are not perfectly nested")
-
- ## aggregate response matrix
- fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
- TRUE else FALSE
-# if (!fullgamma && !global)
-# warning("gamma diversity value might be meaningless")
- ftmp <- vector("list", nlevs)
- for (i in 1:nlevs) {
- ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
- }
-
- ## is there a method/burnin/thin in ... ?
- method <- if (is.null(list(...)$method))
- "r2dtable" else list(...)$method
- burnin <- if (is.null(list(...)$burnin))
- 0 else list(...)$burnin
- thin <- if (is.null(list(...)$thin))
- 1 else list(...)$thin
-
- ## evaluate other arguments
- index <- match.arg(index)
- divfun <- switch(index,
- "renyi" = function(x) renyi(x, scales=scales, hill = TRUE),
- "tsallis" = function(x) tsallis(x, scales=scales, hill = TRUE))
-
- ## cluster membership determination
- nrhs <- rhs
- nrhs <- sapply(nrhs, as.numeric)
- idcl <- function(i) {
- h <- nrhs[,i]
- l <- nrhs[,(i-1)]
- sapply(unique(l), function(i) h[l==i][1])
- }
- id <- lapply(2:nlevs, idcl)
-
- ## this is the function passed to oecosimu
- if (global) {
- wdivfun <- function(x) {
- if (fullgamma) {
- tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- tmp[[nlevs]] <- matrix(colSums(lhs), nrow = 1, ncol = ncol(lhs))
- } else {
- tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- }
- 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 {
- wdivfun <- function(x) {
- if (fullgamma) {
- tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- tmp[[nlevs]] <- matrix(colSums(lhs), nrow = 1, ncol = ncol(lhs))
- } else {
- tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
- }
- a <- sapply(1:nlevs, function(i) divfun(tmp[[i]]))
- am <- lapply(1:(nlevs-1), function(i) {
- sapply(1:length(unique(id[[i]])), function(ii) {
- mean(a[[i]][id[[i]]==ii])
- })
- })
- 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))
- }
- }
- if (nsimul > 0) {
- sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
- burnin=burnin, thin=thin)
- } 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
- attr(sim, "call") <- match.call()
- attr(sim, "index") <- index
- attr(sim, "scales") <- scales
- attr(sim, "global") <- TRUE
- attr(sim, "n.levels") <- nlevs
- attr(sim, "terms") <- tlab
- attr(sim, "model") <- rhs
- class(sim) <- c("multipart", "list")
- sim
+ UseMethod("multipart")
}
Copied: branches/2.0/R/multipart.default.R (from rev 2143, pkg/vegan/R/multipart.default.R)
===================================================================
--- branches/2.0/R/multipart.default.R (rev 0)
+++ branches/2.0/R/multipart.default.R 2012-05-06 13:55:15 UTC (rev 2151)
@@ -0,0 +1,135 @@
+multipart.default <-
+function(y, x, index=c("renyi", "tsallis"), scales = 1,
+ global = FALSE, relative = FALSE, nsimul=99, ...)
+{
+ if (length(scales) > 1)
+ stop("length of 'scales' must be 1")
+ ## evaluate formula
+ lhs <- as.matrix(y)
+ if (missing(x))
+ x <- cbind(level_1=seq_len(nrow(lhs)),
+ leve_2=rep(1, nrow(lhs)))
+ rhs <- data.frame(x)
+ rhs[] <- lapply(rhs, as.factor)
+ rhs[] <- lapply(rhs, droplevels)
+ nlevs <- ncol(rhs)
+ if (nlevs < 2)
+ stop("provide at least two level hierarchy")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vegan -r 2151
More information about the Vegan-commits
mailing list