[Vegan-commits] r2274 - in branches/2.0: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 10 10:00:58 CEST 2012
Author: jarioksa
Date: 2012-09-10 10:00:57 +0200 (Mon, 10 Sep 2012)
New Revision: 2274
Added:
branches/2.0/R/hierParseFormula.R
Removed:
branches/2.0/R/print.adipart.R
branches/2.0/R/print.hiersimu.R
branches/2.0/R/print.multipart.R
Modified:
branches/2.0/NAMESPACE
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
branches/2.0/R/oecosimu.R
branches/2.0/R/print.oecosimu.R
branches/2.0/inst/ChangeLog
branches/2.0/man/adipart.Rd
branches/2.0/man/multipart.Rd
branches/2.0/man/vegan-internal.Rd
Log:
merge r2227:2235 adipart, hiersimu, multipart refactoring (formula, oecosimu)
Modified: branches/2.0/NAMESPACE
===================================================================
--- branches/2.0/NAMESPACE 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/NAMESPACE 2012-09-10 08:00:57 UTC (rev 2274)
@@ -276,7 +276,6 @@
# print: base
S3method(print, CCorA)
S3method(print, MOStest)
-S3method(print, adipart)
S3method(print, adonis)
S3method(print, anosim)
S3method(print, betadisper)
@@ -288,7 +287,6 @@
S3method(print, envfit)
S3method(print, factorfit)
S3method(print, fisherfit)
-S3method(print, hiersimu)
S3method(print, humpfit)
S3method(print, isomap)
S3method(print, mantel)
@@ -297,7 +295,6 @@
S3method(print, monoMDS)
S3method(print, mrpp)
S3method(print, mso)
-S3method(print, multipart)
S3method(print, nestedchecker)
S3method(print, nesteddisc)
S3method(print, nestedn0)
Modified: branches/2.0/R/adipart.default.R
===================================================================
--- branches/2.0/R/adipart.default.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/adipart.default.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -21,7 +21,7 @@
colnames(rhs) <- paste("level", 1:nlevs, sep="_")
tlab <- colnames(rhs)
- ## part check proper design of the model frame
+ ## 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")
@@ -98,12 +98,14 @@
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
+ call <- match.call()
+ call[[1]] <- as.name("adipart")
+ attr(sim, "call") <- call
+ attr(sim$oecosimu$simulated, "index") <- index
+ attr(sim$oecosimu$simulated, "weights") <- weights
attr(sim, "n.levels") <- nlevs
attr(sim, "terms") <- tlab
attr(sim, "model") <- rhs
- class(sim) <- c("adipart", "list")
+ class(sim) <- c("adipart", class(sim))
sim
}
Modified: branches/2.0/R/adipart.formula.R
===================================================================
--- branches/2.0/R/adipart.formula.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/adipart.formula.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,113 +1,19 @@
-adipart.formula <-
-function(formula, data, index=c("richness", "shannon", "simpson"),
- weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
+`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")
+ tmp <- hierParseFormula(formula, data)
+ lhs <- tmp$lhs
+ rhs <- tmp$rhs
- ## 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")
+ ## run simulations
+ sim <- adipart.default(lhs, rhs, index = index, weights = weights,
+ relative = relative, nsimul = nsimul, ...)
+ call <- match.call()
+ call[[1]] <- as.name("adipart")
+ attr(sim, "call") <- call
sim
}
Copied: branches/2.0/R/hierParseFormula.R (from rev 2232, pkg/vegan/R/hierParseFormula.R)
===================================================================
--- branches/2.0/R/hierParseFormula.R (rev 0)
+++ branches/2.0/R/hierParseFormula.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -0,0 +1,20 @@
+"hierParseFormula" <-
+function (formula, data)
+{
+ lhs <- formula[[2]]
+ if (any(attr(terms(formula, data = data), "order") > 1))
+ stop("interactions are not allowed")
+ lhs <- as.matrix(eval(lhs, data))
+ formula[[2]] <- NULL
+ rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
+ rhs[] <- lapply(rhs, function(u) {
+ if (!is.factor(u))
+ u <- factor(u)
+ u
+ })
+ if (length(rhs) < 2)
+ stop("at least 2 hierarchy levels are needed")
+ attr(rhs, "terms") <- NULL
+ list(lhs=lhs, rhs=rhs)
+}
+
Modified: branches/2.0/R/hiersimu.default.R
===================================================================
--- branches/2.0/R/hiersimu.default.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/hiersimu.default.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -15,7 +15,7 @@
colnames(rhs) <- paste("level", 1:nlevs, sep="_")
tlab <- colnames(rhs)
- ## part check proper design of the model frame
+ ## 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")
@@ -82,12 +82,14 @@
# 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()
+ call <- match.call()
+ call[[1]] <- as.name("hiersimu")
+ attr(sim, "call") <- 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")
+ class(sim) <- c("hiersimu", class(sim))
sim
}
Modified: branches/2.0/R/hiersimu.formula.R
===================================================================
--- branches/2.0/R/hiersimu.formula.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/hiersimu.formula.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,96 +1,21 @@
-hiersimu.formula <-
-function(formula, data, FUN, location = c("mean", "median"),
-relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
+`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)
+ tmp <- hierParseFormula(formula, data)
+ lhs <- tmp$lhs
+ rhs <- tmp$rhs
- ## 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(x), nrow = 1, ncol = ncol(x))
- } 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")
+ ## run simulations
+ sim <- hiersimu.default(lhs, rhs, FUN = FUN, location = location,
+ relative = relative, drop.highest = drop.highest,
+ nsimul = nsimul, ...)
+ call <- match.call()
+ call[[1]] <- as.name("hiersimu")
+ attr(sim, "call") <- call
sim
}
+
Modified: branches/2.0/R/multipart.default.R
===================================================================
--- branches/2.0/R/multipart.default.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/multipart.default.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,6 +1,6 @@
-multipart.default <-
-function(y, x, index=c("renyi", "tsallis"), scales = 1,
- global = FALSE, relative = FALSE, nsimul=99, ...)
+`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")
@@ -23,7 +23,7 @@
colnames(rhs) <- paste("level", 1:nlevs, sep="_")
tlab <- colnames(rhs)
- ## part check proper design of the model frame
+ ## 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")
@@ -123,13 +123,15 @@
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
+ call <- match.call()
+ call[[1]] <- as.name("multipart")
+ attr(sim, "call") <- call
+ attr(sim$oecosimu$simulated, "index") <- index
+ attr(sim$oecosimu$simulated, "scales") <- scales
+ attr(sim$oecosimu$simulated, "global") <- TRUE
attr(sim, "n.levels") <- nlevs
attr(sim, "terms") <- tlab
attr(sim, "model") <- rhs
- class(sim) <- c("multipart", "list")
+ class(sim) <- c("multipart", class(sim))
sim
}
Modified: branches/2.0/R/multipart.formula.R
===================================================================
--- branches/2.0/R/multipart.formula.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/multipart.formula.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,139 +1,20 @@
-multipart.formula <-
-function(formula, data, index=c("renyi", "tsallis"), scales = 1,
- global = FALSE, relative = FALSE, nsimul=99, ...)
+`multipart.formula` <-
+ function(formula, data, 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 <- 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")
+ tmp <- hierParseFormula(formula, data)
+ lhs <- tmp$lhs
+ rhs <- tmp$rhs
- ## 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(x), nrow = 1, ncol = ncol(x))
- } 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(x), nrow = 1, ncol = ncol(x))
- } 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")
+ ## run simulations
+ sim <- multipart.default(lhs, rhs, index = index, scales = scales,
+ global = global, relative = relative,
+ nsimul = nsimul, ...)
+ call <- match.call()
+ call[[1]] <- as.name("multipart")
+ attr(sim, "call") <- call
sim
}
Modified: branches/2.0/R/oecosimu.R
===================================================================
--- branches/2.0/R/oecosimu.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/oecosimu.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -128,6 +128,7 @@
ind$oecosimu <- list(z = z, means = means, pval = p, simulated=simind,
method=method,
statistic = indstat, alternative = alternative)
+ attr(ind, "call") <- match.call()
class(ind) <- c("oecosimu", class(ind))
ind
}
Deleted: branches/2.0/R/print.adipart.R
===================================================================
--- branches/2.0/R/print.adipart.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/print.adipart.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,32 +0,0 @@
-print.adipart <-
-function(x, ...)
-{
- n <- if (is.null(x$oecosimu$simulated))
- 0 else ncol(x$oecosimu$simulated)
- if (n > 0)
- cat("adipart with", n, "simulations using method",
- dQuote(x$oecosimu$method), "\n")
- else
- cat("adipart ")
- att <- attributes(x)
- att$names <- att$call <- att$class <- att$n.levels <- att$terms <- att$model <- NULL
- cat("with", paste(names(att), att, collapse=", "))
-
- cat("\n\n")
- cl <- class(x)
- if (length(cl) > 1 && cl[2] != "list") {
- NextMethod("print", x)
- cat("\n")
- }
- 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)
- printCoefmat(m, ...)
- invisible(x)
-}
Deleted: branches/2.0/R/print.hiersimu.R
===================================================================
--- branches/2.0/R/print.hiersimu.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/print.hiersimu.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,16 +0,0 @@
-print.hiersimu <-
-function (x, ...)
-{
- cat("hiersimu with", ncol(x$oecosimu$simulated), "simulations\n\n")
- cl <- class(x)
- if (length(cl) > 1 && cl[2] != "list") {
- NextMethod("print", x)
- cat("\n")
- }
- qu <- apply(x$oecosimu$simulated, 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)
- printCoefmat(m, ...)
- invisible(x)
-}
Deleted: branches/2.0/R/print.multipart.R
===================================================================
--- branches/2.0/R/print.multipart.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/print.multipart.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,27 +0,0 @@
-print.multipart <-
-function(x, ...)
-{
- 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=", "))
- cat("\n\n")
- cl <- class(x)
- if (length(cl) > 1 && cl[2] != "list") {
- NextMethod("print", x)
- cat("\n")
- }
- 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)
- printCoefmat(m, ...)
- invisible(x)
-}
Modified: branches/2.0/R/print.oecosimu.R
===================================================================
--- branches/2.0/R/print.oecosimu.R 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/R/print.oecosimu.R 2012-09-10 08:00:57 UTC (rev 2274)
@@ -1,12 +1,16 @@
`print.oecosimu` <-
function(x, ...)
{
+ xx <- x ## return unmodified input object
attr(x$oecosimu$method, "permfun") <- NULL
- cat("oecosimu with", ncol(x$oecosimu$simulated), "simulations\n")
- cat("simulation method", x$oecosimu$method)
+ cat(as.character(attr(x,"call")[[1]]), "object\n\n")
+ writeLines(strwrap(pasteCall(attr(x, "call"))))
+ cat("\n")
+ cat("simulation method", x$oecosimu$method, "with",
+ ncol(x$oecosimu$simulated), "simulations\n")
if (length(att <- attributes(x$oecosimu$simulated)) > 1) {
att$dim <- NULL
- cat(" with", paste(names(att), att, collapse=", "))
+ cat("options: ", paste(names(att), att, collapse=", "))
}
alt.char <- switch(x$oecosimu$alternative,
two.sided = "not equal to",
@@ -17,9 +21,10 @@
cat("\n\n")
cl <- class(x)
- if (length(cl) > 1 && cl[2] != "list") {
- NextMethod("print", x)
- cat("\n")
+ if ((length(cl) > 1 && cl[2] != "list" ) &&
+ !any(cl %in% c("adipart", "hiersimu", "multipart"))) {
+ NextMethod("print", x)
+ cat("\n")
}
probs <- switch(x$oecosimu$alternative,
two.sided = c(0.025, 0.5, 0.975),
@@ -35,5 +40,7 @@
cat("\nNumber of NA cases removed from simulations:\n",
nacount, "\n")
}
- invisible(x)
+ invisible(xx)
}
+
+
Modified: branches/2.0/inst/ChangeLog
===================================================================
--- branches/2.0/inst/ChangeLog 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/inst/ChangeLog 2012-09-10 08:00:57 UTC (rev 2274)
@@ -14,6 +14,11 @@
* merge r2244: more portable doc/Makefile
* merge r2237 thru 2240: add labels.envfit() and "labels" arg to
plot.envfit().
+ * merge r2227 thru 2235: adipart, hiersimu, multipart code
+ refactoring (especially formula method) and making to inherit from
+ oecosimu in printing the results. The merge did not apply quite
+ cleanly, but oecosimu, print.oecosimu and NAMESPACE needed manual
+ editing (beware).
* merge r2225: biplot.rda 'type' fix.
Version 2.0-4 (released June 18, 2012)
Modified: branches/2.0/man/adipart.Rd
===================================================================
--- branches/2.0/man/adipart.Rd 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/man/adipart.Rd 2012-09-10 08:00:57 UTC (rev 2274)
@@ -3,11 +3,10 @@
\alias{adipart}
\alias{adipart.default}
\alias{adipart.formula}
-\alias{print.adipart}
\alias{hiersimu}
\alias{hiersimu.default}
\alias{hiersimu.formula}
-\alias{print.hiersimu}
+
\title{Additive Diversity Partitioning and Hierarchical Null Model Testing}
\description{
In additive diversity partitioning, mean values of alpha diversity at lower levels of a sampling
@@ -37,10 +36,9 @@
all rows are in the same group in the second level.}
\item{formula}{A two sided model formula in the form \code{y ~ x}, where \code{y}
is the community data matrix with samples as rows and species as column. Right
- hand side (\code{x}) must contain factors referring to levels of sampling hierarchy,
+ hand side (\code{x}) must grouping vaiables referring to levels of sampling hierarchy,
terms from right to left will be treated as nested (first column is the lowest,
- last is the highest level). These variables must be factors in order to unambiguous
- handling. Interaction terms are not allowed.}
+ last is the highest level, at least two levels specified). Interaction terms are not allowed.}
\item{data}{A data frame where to look for variables defined in the right hand side
of \code{formula}. If missing, variables are looked in the global environment.}
\item{index}{Character, the diversity index to be calculated (see Details).}
@@ -128,10 +126,10 @@
out <- rep(1, length(x))
for (i in 2:(length(cut) - 1))
out[which(x > cut[i] & x <= cut[(i + 1)])] <- i
- return(as.factor(out))}
+ return(out)}
## The hierarchy of sample aggregation
levsm <- data.frame(
- l1=as.factor(1:nrow(mite)),
+ l1=1:nrow(mite),
l2=cutter(mite.xy$y, cut = seq(0, 10, by = 2.5)),
l3=cutter(mite.xy$y, cut = seq(0, 10, by = 5)),
l4=cutter(mite.xy$y, cut = seq(0, 10, by = 10)))
Modified: branches/2.0/man/multipart.Rd
===================================================================
--- branches/2.0/man/multipart.Rd 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/man/multipart.Rd 2012-09-10 08:00:57 UTC (rev 2274)
@@ -3,7 +3,6 @@
\alias{multipart}
\alias{multipart.default}
\alias{multipart.formula}
-\alias{print.multipart}
\title{Multiplicative Diversity Partitioning}
\description{
@@ -26,10 +25,9 @@
all rows are in the same group in the second level.}
\item{formula}{A two sided model formula in the form \code{y ~ x}, where \code{y}
is the community data matrix with samples as rows and species as column. Right
- hand side (\code{x}) must contain factors referring to levels of sampling hierarchy,
+ hand side (\code{x}) must grouping vaiables referring to levels of sampling hierarchy,
terms from right to left will be treated as nested (first column is the lowest,
- last is the highest level). These variables must be factors in order to unambiguous
- handling. Interaction terms are not allowed.}
+ last is the highest level, at least two levels specified). Interaction terms are not allowed.}
\item{data}{A data frame where to look for variables defined in the right hand side
of \code{formula}. If missing, variables are looked in the global environment.}
\item{index}{Character, the entropy index to be calculated (see Details).}
@@ -105,10 +103,10 @@
out <- rep(1, length(x))
for (i in 2:(length(cut) - 1))
out[which(x > cut[i] & x <= cut[(i + 1)])] <- i
- return(as.factor(out))}
+ return(out)}
## The hierarchy of sample aggregation
levsm <- data.frame(
- l1=as.factor(1:nrow(mite)),
+ l1=1:nrow(mite),
l2=cutter(mite.xy$y, cut = seq(0, 10, by = 2.5)),
l3=cutter(mite.xy$y, cut = seq(0, 10, by = 5)),
l4=cutter(mite.xy$y, cut = seq(0, 10, by = 10)))
Modified: branches/2.0/man/vegan-internal.Rd
===================================================================
--- branches/2.0/man/vegan-internal.Rd 2012-09-10 07:02:23 UTC (rev 2273)
+++ branches/2.0/man/vegan-internal.Rd 2012-09-10 08:00:57 UTC (rev 2274)
@@ -10,6 +10,7 @@
\alias{ordiArrowMul}
\alias{ordiArgAbsorber}
\alias{veganCovEllipse}
+\alias{hierParseFormula}
\title{Internal vegan functions}
@@ -31,6 +32,7 @@
permuted.index(n, strata)
pasteCall(call, prefix = "Call:")
veganCovEllipse(cov, center = c(0, 0), scale = 1, npoints = 100)
+hierParseFormula(formula, data)
}
\details{ The description here is only intended for \pkg{vegan}
@@ -79,6 +81,11 @@
\code{veganCovEllipse} finds the coordinates for drawing a
covariance ellipse.
+
+ \code{hierParseFormula} returns a list of one matrix (left hand side)
+ and a model frame with factors representing hierarchy levels
+ (right hand side) to be used in \code{\link{adipart}},
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vegan -r 2274
More information about the Vegan-commits
mailing list