[Vegan-commits] r2714 - in pkg/vegan: R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 20 12:26:37 CET 2013
Author: jarioksa
Date: 2013-11-20 12:26:37 +0100 (Wed, 20 Nov 2013)
New Revision: 2714
Added:
pkg/vegan/R/anovacca.byterm.R
Modified:
pkg/vegan/R/anovacca.R
pkg/vegan/inst/ChangeLog
Log:
add by='term' and by='margin' cases to new anova(.)cca tests
Modified: pkg/vegan/R/anovacca.R
===================================================================
--- pkg/vegan/R/anovacca.R 2013-11-19 14:40:44 UTC (rev 2713)
+++ pkg/vegan/R/anovacca.R 2013-11-20 11:26:37 UTC (rev 2714)
@@ -30,19 +30,19 @@
## stop permutations block
## see if this was a list of ordination objects
dotargs <- list(...)
+ ## we do not want to give dotargs to anova.ccalist, but we
+ ## evaluate 'parallel' and 'model' here
+ if (is.null(dotargs$model))
+ model <- "reduced"
+ else
+ model <- dotargs$model
+ if (is.null(dotargs$parallel))
+ parallel <- NULL
+ else
+ parallel <- dotargs$parallel
if (length(dotargs)) {
isCCA <- sapply(dotargs, function(z) inherits(z, "cca"))
if (any(isCCA)) {
- ## we do not want to give dotargs to anova.ccalist, but we
- ## evaluate 'parallel' and 'model' here
- if (is.null(dotargs$model))
- model <- "reduced"
- else
- model <- dotargs$model
- if (is.null(dotargs$parallel))
- parallel <- NULL
- else
- parallel <- dotargs$parallel
dotargs <- dotargs[isCCA]
object <- c(list(object), dotargs)
sol <-
@@ -55,8 +55,15 @@
}
## by cases
if (!is.null(by)) {
- by <- match.arg(by, c("axis", "terms", "margin"))
- .NotYetUsed("by")
+ by <- match.arg(by, c("terms", "margin"))
+ sol <- switch(by,
+ "terms" = anovacca.byterm(object,
+ permutations = permutations,
+ model = model, parallel = parallel),
+ "margin" = anovacca.bymargin(object,
+ permutations = permutations,
+ model = model, parallel = parallel))
+ return(sol)
}
## basic overall test
tst <- permutest.cca(object, permutations = permutations, ...)
Added: pkg/vegan/R/anovacca.byterm.R
===================================================================
--- pkg/vegan/R/anovacca.byterm.R (rev 0)
+++ pkg/vegan/R/anovacca.byterm.R 2013-11-20 11:26:37 UTC (rev 2714)
@@ -0,0 +1,69 @@
+### Implementation of by-cases for vegan 2.2 versions of
+### anova.cca. These are all internal functions that are not intended
+### to be called by users in normal sessions, but they should be
+### called from anova.cca (2.2). Therefore the user interface is rigid
+### and input is not checked. The 'permutations' should be a
+### permutation matrix.
+
+### by = terms builds models as a sequence of adding terms and submits
+### this to anova.ccalist
+
+`anovacca.byterm` <-
+ function(object, permutations, model, parallel)
+{
+ ## We need term labels but without Condition() terms
+ trms <- terms(object)
+ trmlab <- attr(trms, "term.labels")
+ trmlab <- trmlab[trmlab %in% attr(terms(object$terminfo),
+ "term.labels")]
+ m0 <- update(object, . ~ 1)
+ mods <- list(m0)
+ for(i in seq_along(trmlab)) {
+ fla <- paste(". ~ . + ", trmlab[i])
+ mods[[i+1]] <- update(mods[[i]], fla)
+ }
+ ## The result. Should be reformatted
+ anova.ccalist(mods, permutations = permutations, model = model, parallel = parallel)
+}
+
+## by = margin: this is not a anova.ccalist case, but we omit each
+## term in turn and compare against the complete model.
+
+`anovacca.bymargin` <-
+ function(object, permutations, ...)
+{
+ nperm <- nrow(permutations)
+ ## We need term labels but without Condition() terms
+ trms <- terms(object)
+ trmlab <- attr(trms, "term.labels")
+ trmlab <- trmlab[trmlab %in% attr(terms(object$terminfo),
+ "term.labels")]
+ ## baseline: all terms
+ big <- permutest(object, permutations, ...)
+ dfbig <- big$df[2]
+ chibig <- big$chi[2]
+ scale <- big$den/dfbig
+ ## Collect all marginal models
+ mods <- lapply(trmlab, function(nm, ...)
+ permutest(update(object, paste(".~.-", nm)),
+ permutations, ...))
+ ## Chande in df
+ Df <- sapply(mods, function(x) x$df[2]) - dfbig
+ ## F of change
+ Chisq <- sapply(mods, function(x) x$chi[2]) - chibig
+ Fstat <- (Chisq/Df)/(chibig/dfbig)
+ ## Simulated F-values
+ Fval <- sapply(mods, function(x) x$den)
+ Fval <- sweep(Fval, 1, big$den, "-")
+ Fval <- sweep(Fval, 2, Df, "/")
+ Fval <- sweep(Fval, 1, scale, "/")
+ ## Simulated P-values
+ Pval <- (colSums(sweep(Fval, 2, Fstat, ">=")) + 1)/(nperm + 1)
+ ## Collect results to anova data.frame
+ out <- data.frame(c(Df, dfbig), c(Chisq, chibig),
+ c(Fstat, NA), c(Pval, NA))
+ colnames(out) <- c("Df", "Chisq", "F", "Pr(>F)")
+ rownames(out) <- c(trmlab, "Residual")
+ class(out) <- c("anova", "data.frame")
+ out
+}
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2013-11-19 14:40:44 UTC (rev 2713)
+++ pkg/vegan/inst/ChangeLog 2013-11-20 11:26:37 UTC (rev 2714)
@@ -10,9 +10,16 @@
the function no more adapts the number of iterations for the
P-value, and arguments 'step' and 'perm.max' were removed.
Instead, permute package is used to create a permutation matrix
- used in all cases with fixed number of permutations. At the first
- stage, only the overall test is provided.
+ used in all cases with fixed number of permutations.
+ In addition to the overall test, the function allows now testing
+ a sequence of models (anova.ccalist). Specific tests provided are
+ by = "term" (anovacca.byterm) which is fitted as a sequence of
+ models in anova.ccalist. Case by = "margin" directly calls
+ permutest.cca and gets the significances from differences of
+ residual variation similarly as anova.ccalist. Case by = "axis" is
+ not yet implemented.
+
* commsim: documentation (commsim.Rd) was restructured so that
nullmodels were collected under separate sections with a brief
introductory text and shorter specific text of the
More information about the Vegan-commits
mailing list