[Vegan-commits] r1273 - in pkg/vegan: R inst man

[noreply at r-forge.r-project.org]

Author: jarioksa
Date: 2010-08-20 14:47:55 +0200 (Fri, 20 Aug 2010)
New Revision: 1273

Added:
   pkg/vegan/R/anova.ccalist.R
Modified:
   pkg/vegan/inst/ChangeLog
   pkg/vegan/man/anova.cca.Rd
Log:
fist try on anova.ccalist or comparison of several cca/rda/capscace models

Added: pkg/vegan/R/anova.ccalist.R
===================================================================
--- pkg/vegan/R/anova.ccalist.R	                        (rev 0)
+++ pkg/vegan/R/anova.ccalist.R	2010-08-20 12:47:55 UTC (rev 1273)
@@ -0,0 +1,45 @@
+### A test of concept function for significance test of two (or more)
+### cca, rda or capscale models. A "test of concept" means that there
+### are no sanity chekcs, but only the basic code with assumption that
+### input is correct.
+
+### The idea is to run permutest with the same random number seed and
+### take the difference of permutation results (denominator and
+### numerator of pseudo-F) and take P-values of that.  Assume the two
+### models are m1 <- cca(y ~ x1) and m2 <- cca(y ~ x1 + x2). Currently
+### the difference can be analysed as a partial model, where the
+### common part is partialled out m12 <- cca(y ~ Condition(x1) +
+### x2). The direct difference of two models suggested here does not
+### produce the same results except with model = "direct", because the
+### partial model will permute residuals after conditions.
+`anova.ccalist` <- function(object, ...)
+{
+    objects <- list(object, ...)
+    nmodels <- length(objects)
+    ## Collect statistics
+    resdf <- sapply(objects, function(x) x$CA$rank)
+    resdev <- sapply(objects, function(x) x$CA$tot.chi)
+    ## residual df (df) do not necessarily decrease with rank deficit
+    ## models: we also need model (CCA) df and their change
+    moddf <- sapply(objects, function(x) x$CCA$qrank)
+    ## ANOVA table except test
+    table <- data.frame(moddf, resdf, resdev, c(NA, diff(moddf)),
+                        c(NA, -diff(resdev)))
+    dimnames(table) <- list(1:nmodels, c("Model.Df", "Res.Df", "RSS", "Df", "Sum of Sq"))
+    variables <- sapply(objects,
+                        function(x) paste(deparse(formula(x)), collapse="\n"))
+    title <- "Analysis of Variance Table\n"
+    topnote <- paste("Model ", format(1:nmodels), ": ", variables,
+                     sep = "", collapse = "\n")
+    table <- structure(table, heading = c(title, topnote))
+    ##Collect tests
+    mods <- list()
+    for(i in 1:nmodels) {
+        if (i > 1)
+            assign(".Random.seed", mods[[1]]$Random.seed,
+                   envir = .GlobalEnv)
+        mods[[i]] <- permutest(objects[[i]])
+    }
+    class(table) <- c("anova", class(table))
+    table
+}

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2010-08-20 08:08:41 UTC (rev 1272)
+++ pkg/vegan/inst/ChangeLog	2010-08-20 12:47:55 UTC (rev 1273)
@@ -7,8 +7,29 @@
 	* New version opened with the release of vegan_1.17-4 on August
 	20, 2010.
 
-	* permustest.cca: defaults to 99 permutations instead of 100 (to
+	* permutest.cca: defaults to 99 permutations instead of 100 (to
 	be consistent).
+
+	* anova.ccalist: Started to work on the concept of permutation
+	test comparison of several cca/ rda/ capscale models. The first
+	stage of the function was added to the R-Forge (undocumented!).
+	Currently the function only collects the statistics to an ANOVA
+	table but performs no tests. My first idea is to run permutest.cca
+	with the same RNG seed and collect the permutations from the
+	differences of these models. I know now that this won't work
+	because permuted residuals differ from the corresponding partial
+	model (except with model = "direct").  Another option that I have
+	on my mind is that I use model.matrix.cca to reconstruct a matrix
+	of RHS of each model, and then run a series partial model
+	ANOVAs. Here the problem is that the statistics will be
+	sequential, i.e., the residual of each model is from the next
+	model and not from the final model (bigmodel). Naturally, we have
+	to figure out how the permutation should be done. Probably a note
+	to R-Devel to open up discussion on the issue. NOTE: there is no
+	automatic delegation to this function. If you try this, you must
+	explcititly write "anova.ccalist(model1, model2, model3). NOTE2:
+	There are no tests of the sanity of the input: be careful. NOTE3:
+	No test are provided: still under construction.
 	
 Version 1.18-9 (closed August 20, 2010)
 

Modified: pkg/vegan/man/anova.cca.Rd
===================================================================
--- pkg/vegan/man/anova.cca.Rd	2010-08-20 08:08:41 UTC (rev 1272)
+++ pkg/vegan/man/anova.cca.Rd	2010-08-20 12:47:55 UTC (rev 1273)
@@ -8,6 +8,7 @@
 \alias{permutest.default}
 \alias{permutest.cca}
 \alias{print.permutest.cca}
+\alias{anova.ccalist}
 
 \title{Permutation Test for Constrained Correspondence Analysis,
   Redundancy Analysis and Constrained Analysis of Principal Coordinates }