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

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

Author: jarioksa
Date: 2009-10-05 21:39:01 +0200 (Mon, 05 Oct 2009)
New Revision: 1032

Added:
   pkg/vegan/R/ordistep.R
Modified:
   pkg/vegan/inst/ChangeLog
   pkg/vegan/man/add1.cca.Rd
Log:
ordistep: new function to model choice in the cca family based on permutation P values (alternative to step using AIC)

Added: pkg/vegan/R/ordistep.R
===================================================================
--- pkg/vegan/R/ordistep.R	                        (rev 0)
+++ pkg/vegan/R/ordistep.R	2009-10-05 19:39:01 UTC (rev 1032)
@@ -0,0 +1,60 @@
+`ordistep` <-
+    function(object, scope, Pin = 0.05, Pout = 0.1,
+             pstep = 100, perm.max = 1000, steps=50, trace = TRUE, ...)
+{
+    if (!inherits(object, "cca"))
+        stop("function can be only used to 'cca' and related objects")
+    if (missing(scope))
+        stop("absolutely needs scope")
+    mod <- eval.parent(update(object, . ~  1))
+    for (i in 1:steps){
+        change <- NULL
+        ## Consider dropping
+        if (length(drop.scope(mod))) {
+            aod <- drop1(mod, test="perm", pstep = pstep, perm.max = perm.max,
+                         alpha = Pout)
+            aod <- aod[-1,]
+            o <- order(-aod[,5], aod[,4], aod[,2])
+            aod <- aod[o,]
+            rownames(aod) <- paste("-", rownames(aod), sep = " ")
+            if (trace) {
+                cat("\n")
+                print(aod)
+            }
+            if (aod[1,5] > Pout) {
+                change <- rownames(aod)[1]
+                mod <- eval.parent(update(mod, paste("~  .", change)))
+                if (trace) {
+                    cat("\n")
+                    print(mod$call)
+                }
+            }
+        }
+        ## Consider adding
+        if (length(add.scope(formula(mod), scope))) {
+            aod <- add1(mod, scope = scope, test = "perm", pstep = pstep,
+                        perm.max = perm.max, alpha = Pin)
+            aod <- aod[-1,]
+            o <- order(aod[,5], aod[,4], aod[,2])
+            aod <- aod[o,]
+            rownames(aod) <- paste("+", rownames(aod), sep = " ")
+            if (trace) {
+                cat("\n")
+                print(aod)
+            }
+            if (aod[1,5] <= Pin) {
+                change <- rownames(aod)[1]
+                mod <- eval.parent(update(mod, paste( "~  .",change)))
+                if (trace) {
+                    cat("\n")
+                    print(mod$call)
+                }
+            }
+        }
+        ## No drop, no add: done
+        if (is.null(change))
+            break
+    }
+    cat("\n")
+    mod
+}

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2009-10-05 05:09:17 UTC (rev 1031)
+++ pkg/vegan/inst/ChangeLog	2009-10-05 19:39:01 UTC (rev 1032)
@@ -4,6 +4,17 @@
 
 Version 1.16-31 (opened October 3, 2009)
 
+	* ordistep: new function for automatic model choice in the cca
+	family of methods. This is similar to step(), but uses permutation
+	P-values instead of AIC (that does not exist for these methods,
+	and sometimes fails). This is similar to forward.sel() in Stephane
+	Dray's 'packfor' package at R-Forge. Better things here are that
+	the function has both drop and add instead of only forward
+	selection, handles correctly factors, and works with cca(), rda()
+	and capscale(). On the other hand, it is much slower than
+	forward.sel(). The 'scope' needs still work. Not yet documented
+	(alias with add1.cca). 
+
 	* ordisurf: pushes out limits of convex hull so that contours
 	better cover corners and margins of points.
 

Modified: pkg/vegan/man/add1.cca.Rd
===================================================================
--- pkg/vegan/man/add1.cca.Rd	2009-10-05 05:09:17 UTC (rev 1031)
+++ pkg/vegan/man/add1.cca.Rd	2009-10-05 19:39:01 UTC (rev 1032)
@@ -1,6 +1,7 @@
 \name{add1.cca}
 \alias{add1.cca}
 \alias{drop1.cca}
+\alias{ordistep}
 
 \title{Add or Drop  Single Terms to a Constrained Ordination Model }
 \description{