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

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

Author: jarioksa
Date: 2010-08-23 16:46:53 +0200 (Mon, 23 Aug 2010)
New Revision: 1274

Modified:
   pkg/vegan/R/anova.ccalist.R
   pkg/vegan/inst/ChangeLog
Log:
more or less working version of anova.ccalist

Modified: pkg/vegan/R/anova.ccalist.R
===================================================================
--- pkg/vegan/R/anova.ccalist.R	2010-08-20 12:47:55 UTC (rev 1273)
+++ pkg/vegan/R/anova.ccalist.R	2010-08-23 14:46:53 UTC (rev 1274)
@@ -17,21 +17,13 @@
     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))
+    N <- nrow(object$CA$u)
+    ranks <- sapply(objects, function(x)
+                    if (is.null(x$CCA$qrank)) 0 else x$CCA$qrank)
+    resdf <- N - ranks - 1
+    resdev <- sapply(objects, deviance)
+    moddev <- c(NA, -diff(resdev))
+    moddf <- c(NA, diff(ranks))
     ##Collect tests
     mods <- list()
     for(i in 1:nmodels) {
@@ -40,6 +32,27 @@
                    envir = .GlobalEnv)
         mods[[i]] <- permutest(objects[[i]])
     }
-    class(table) <- c("anova", class(table))
+    ## Differences of permutations. In permutation F values, numerator
+    ## is taken from each model, but all use the same denominator from
+    ## the largest model.
+    bigmodel <- which.min(resdf)
+    F <- moddev/moddf/resdev[bigmodel]*resdf[bigmodel]
+    den <- mods[[bigmodel]]$den/mods[[bigmodel]]$df[2]
+    Pval <- rep(NA, nmodels)
+    for (i in 2:nmodels) {
+        F.perm <- abs(mods[[i]]$num - mods[[i-1]]$num)/abs(moddf[i])/abs(den)
+        Pval[i] <- (sum(F.perm >= F[i]) + 1)/(length(F.perm) + 1)
+    }   
+    ## ANOVA table 
+    table <- data.frame(resdf, resdev, moddf, moddev, F, Pval)
+    dimnames(table) <- list(1:nmodels, c("Res.Df", "RSS", "Df", "Sum of Sq",
+                                         "F", "Pr(>F)"))
+    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))
+    class(table) <- c("anova.cca", "anova", class(table))
     table
 }

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2010-08-20 12:47:55 UTC (rev 1273)
+++ pkg/vegan/inst/ChangeLog	2010-08-23 14:46:53 UTC (rev 1274)
@@ -13,23 +13,23 @@
 	* 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.
+	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: Does not yet work
+	with Null model (y ~ 1) since permutest.cca won't work with them
+	(anova.ccanull would work).
 	
 Version 1.18-9 (closed August 20, 2010)