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

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

Author: jarioksa
Date: 2013-10-14 09:40:57 +0200 (Mon, 14 Oct 2013)
New Revision: 2633

Added:
   pkg/vegan/R/permutest.ccalist.R
Modified:
   pkg/vegan/DESCRIPTION
   pkg/vegan/inst/ChangeLog
Log:
add permutest.ccalist() to analyse a sequence of cca models (experimental)

Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION	2013-10-03 12:49:48 UTC (rev 2632)
+++ pkg/vegan/DESCRIPTION	2013-10-14 07:40:57 UTC (rev 2633)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 2.1-36
-Date: September 25, 2013
+Version: 2.1-37
+Date: October 14, 2013
 Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, 
    Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, 
    M. Henry H. Stevens, Helene Wagner  

Added: pkg/vegan/R/permutest.ccalist.R
===================================================================
--- pkg/vegan/R/permutest.ccalist.R	                        (rev 0)
+++ pkg/vegan/R/permutest.ccalist.R	2013-10-14 07:40:57 UTC (rev 2633)
@@ -0,0 +1,74 @@
+`permutest.ccalist` <-
+    function(x, ..., permutations = 99)
+{
+    ## Collect cca class objects. FIXME: Eventually this should be in
+    ## a function that calls permutest.ccalist after collecting model
+    ## objects from dotargs.
+    dotargs <- list(...)
+    if (length(dotargs)) {
+        isCCA <- sapply(dotargs, function(z) inherits(z, "cca"))
+        dotargs <- dotargs[isCCA]
+        if (length(dotargs))
+            x <- c(list(x), dotargs)
+    }
+    nmodels <- length(x)
+    ## No. of observations and check
+    N <- sapply(x, nobs)
+    if (!all(N = N[1]))
+        stop("models have different numbers of observations")
+    else
+        N <- N[1]
+    ## Create permutation matrix if it does not exist. FIXME: should
+    ## take arguments for restricted permutation
+    if (length(permutations) == 1) 
+            permutations <- shuffleSet(N, permutations)
+    ## permutations is now matrix
+    nperm <- nrow(permutations)
+    ## check
+    if (ncol(permutations) != N)
+        stop(gettextf("permutation matrix has %d columns, but you have %d sites",
+                      ncol(nperm), N))
+    ## All models are evaluated in permutest.cca with identical
+    ## permutations so that the differences of single permutations can
+    ## be used to assess the significance of differences of fitted
+    ## models. This strictly requires nested models (not checked
+    ## here): all terms of the smaller model must be included in the
+    ## larger model. FIXME: should pass arguments to permutest.cca.
+    mods <- lapply(x, function(z)
+                   permutest(z, permutations = permutations))
+    dfs <- sapply(mods, function(z) z$df)
+    dev <- sapply(mods, function(z) z$chi)
+    resdf <- dfs[2,]
+    df <- -diff(resdf)
+    resdev <- dev[2,]
+    changedev <- -diff(resdev)
+    big <- which.min(resdf)
+    scale <- resdev[big]/resdf[big]
+    fval <- changedev/df/scale
+    ## Collect permutation results: denominator of F varies in each
+    ## permutation.
+    pscale <- mods[[big]]$den/resdf[big]
+    ## Numerator of F
+    pfvals <- sapply(mods, function(z) z$num)
+    pfvals <- apply(pfvals, 1, diff)
+    ## dropped to vector?
+    if (!is.matrix(pfvals))
+        pfvals <- matrix(pfvals, nrow=1, ncol=nperm)
+    pval <- rowSums(sweep(pfvals, 1, fval, ">="))
+    pval <- (pval + 1)/(nperm+1)
+    pfvals <- sweep(pfvals, 1, df, "/")
+    pfvals <- sweep(pfvals, 2, pscale, "/")
+    ## collect table
+    table <- data.frame(resdf, resdev, c(NA, df),
+                        c(NA,changedev), c(NA,fval), c(NA,pval))
+    dimnames(table) <- list(1L:nmodels, c("Resid. Df", "Res. Chisq", 
+                                          "Df", "Chisq", "F", "Pr(>F)"))
+    ## Collect header information
+    formulae <- sapply(x, function(z) deparse(formula(z)))
+    head <- paste0("Permutation tests for ", x[[1]]$method, " under ",
+                  mods[[1]]$model, " model\nwith ", nperm,
+                   " permutations\n")
+    topnote <- paste("Model ", format(1L:nmodels), ": ", formulae,
+                     sep = "", collapse = "\n")
+    structure(table, heading=c(head,topnote), class = c("anova", "data.frame"))
+}

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2013-10-03 12:49:48 UTC (rev 2632)
+++ pkg/vegan/inst/ChangeLog	2013-10-14 07:40:57 UTC (rev 2633)
@@ -2,8 +2,29 @@
 
 VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
 
-Version 2.1-36 (opened September 25, 2013)
+Version 2.1-37 (opened October 14, 2013)
 
+	* permutest.cca: added new function permutest.ccalist() to compare
+	a sequence of models. The function is still experimental ("proof
+	of the concept") and unexported. If this stays in vegan, it should
+	eventually be called from anova.cca() or permutest.cca(). This
+	would bring along a change of API to, say, anova.cca(object, ...,
+	alpha=...): the dots must follow the first argument which turns of
+	positional and partial matching of arguments so that the function
+	can collect the "cca" models. We must decide whether the new
+	function is worth such a change that can make life harder for
+	ordinary users.
+
+	The function is based on calling permutest.cca for each model with
+	identical permutations. We can then compare the change in model
+	for each permutation and collect the test statistics for
+	differences. This requires that the models really are nested so
+	that residual deviance certainly decreases in bigger model
+	(testing theory requires nesting, but this is commonly violated by
+	users: here nesting is necessary). 
+	
+Version 2.1-36 (closed October 14, 2013)
+
 	* opened with the release of vegan 2.0-9.
 
 	* decostand(..., "normalize") uses now .Machine$double.eps to