[Vegan-commits] r1395 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 8 19:10:29 CET 2010
Author: jarioksa
Date: 2010-12-08 19:10:29 +0100 (Wed, 08 Dec 2010)
New Revision: 1395
Added:
branches/1.17/R/ordiR2step.R
Modified:
branches/1.17/inst/ChangeLog
branches/1.17/man/ordistep.Rd
Log:
cp ordiR2step to branches/1.17
Copied: branches/1.17/R/ordiR2step.R (from rev 1394, pkg/vegan/R/ordiR2step.R)
===================================================================
--- branches/1.17/R/ordiR2step.R (rev 0)
+++ branches/1.17/R/ordiR2step.R 2010-12-08 18:10:29 UTC (rev 1395)
@@ -0,0 +1,98 @@
+### Forward selection to maximize R2.adjusted, but stopping once the
+### R2.adjusted of the maximum model ('scope') is exceeded, after
+### Blanchet, Legendre & Borcard: Ecology 89, 2623--2623; 2008.
+
+`ordiR2step` <-
+ function(object, scope, direction = c("both", "forward"),
+ Pin = 0.05, pstep = 100, perm.max = 1000,
+ trace = TRUE, ...)
+{
+ direction <- match.arg(direction)
+ if (missing(scope))
+ stop("needs scope")
+ ## Works only for rda(): cca() does not have (yet) R2.adjusted
+ if (!inherits(object, "rda"))
+ stop("can be used only with rda() or capscale()")
+ ## Get R2 of the original object
+ if (is.null(object$CCA))
+ R2.0 <- 0
+ else
+ R2.0 <- RsquareAdj(object)$adj.r.squared
+ ## Get R2 of the scope
+ if (inherits(scope, "rda"))
+ scope <- delete.response(formula(scope))
+ if (!inherits(scope, "formula"))
+ scope <- reformulate(scope)
+ R2.all <- RsquareAdj(update(object, scope))$adj.r.squared
+ ## Check that the full model can be evaluated
+ if (is.na(R2.all))
+ stop("the upper scope cannot be fitted (too many terms?)")
+ ## Collect data to anotab returned as the 'anova' object
+ anotab <- list()
+ ## Step forward and continue as long as R2.adj improves and R2.adj
+ ## remains below R2.adj < R2.all
+ R2.previous <- R2.0
+ drops <- NULL
+ repeat {
+ if (trace) {
+ cat("Step: R2.adj=", R2.previous, "\n")
+ cat(pasteCall(formula(object)), "\n")
+ }
+ adds <- add.scope(object, scope)
+ if (direction == "both")
+ drops <- drop.scope(object)
+ ## Nothing to add or drop, and we're done: break
+ if (length(adds) == 0 && length(drops) == 0)
+ break
+ R2.adds <- numeric(length(adds) + length(drops))
+ if (length(adds))
+ adds <- paste("+", adds)
+ if (length(drops))
+ drops <- paste("-", drops)
+ names(R2.adds) <- c(adds, drops)
+ ## Loop over add scope
+ for (trm in seq_along(R2.adds)) {
+ fla <- paste(". ~ .", names(R2.adds[trm]))
+ R2.tmp <- RsquareAdj(update(object, fla))$adj.r.squared
+ if (!length(R2.tmp))
+ R2.tmp <- 0
+ R2.adds[trm] <- R2.tmp
+ }
+ best <- which.max(R2.adds)
+ if (trace) {
+ out <- sort(c("<All variables>" = R2.all, "<none>" = R2.previous,
+ R2.adds), decreasing = TRUE)
+ out <- as.matrix(out)
+ colnames(out) <- "R2.adjusted"
+ print(out)
+ cat("\n")
+ }
+ ## See if the best should be kept
+ ## First criterion: R2.adj improves and is still lower or
+ ## equal than for the full model of the scope
+ if (R2.adds[best] > R2.previous && R2.adds[best] <= R2.all) {
+ ## Second criterion: added variable is significant
+ tst <- add1(object, scope = adds[best], test="permu",
+ pstep = pstep, perm.max = perm.max,
+ alpha = Pin, trace = FALSE, ...)
+ if (trace) {
+ print(tst[-1,])
+ cat("\n")
+ }
+ if (tst[,"Pr(>F)"][2] > Pin)
+ break
+ fla <- paste("~ .", adds[best])
+ object <- update(object, fla)
+ R2.previous <- RsquareAdj(object)$adj.r.squared
+ anotab <- rbind(anotab, cbind("R2.adj" = R2.previous, tst[2,]))
+ } else {
+ break
+ }
+ }
+ if (NROW(anotab) > 0) {
+ anotab <- rbind(anotab, "<All variables>" = c(R2.all, rep(NA, 5)))
+ class(anotab) <- c("anova", class(anotab))
+ object$anova <- anotab
+ }
+ object
+}
Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog 2010-12-08 17:55:24 UTC (rev 1394)
+++ branches/1.17/inst/ChangeLog 2010-12-08 18:10:29 UTC (rev 1395)
@@ -4,6 +4,9 @@
Version 1.17-5 (opened November 17, 2010)
+ * copied ordiR2step.R from pkg/vegan and merged docs to
+ ordistep.Rd.
+
* merged r1273, 1275: permutest.cca 99 permutation & saves its
call.
Modified: branches/1.17/man/ordistep.Rd
===================================================================
--- branches/1.17/man/ordistep.Rd 2010-12-08 17:55:24 UTC (rev 1394)
+++ branches/1.17/man/ordistep.Rd 2010-12-08 18:10:29 UTC (rev 1395)
@@ -1,6 +1,7 @@
\name{ordistep}
\Rdversion{1.1}
\alias{ordistep}
+\alias{ordiR2step}
\title{
Choose a Model by Permutation Tests in Constrained Ordination
@@ -8,28 +9,36 @@
\description{
Automatic stepwise model building for constrained ordination methods
(\code{\link{cca}}, \code{\link{rda}}, \code{\link{capscale}}).
- The function is modelled after \code{\link{step}} and can do forward,
- backward and stepwise model selection.
+ The function \code{ordistep} is modelled after \code{\link{step}} and
+ can do forward, backward and stepwise model selection using permutation tests.
+ Function \code{ordiR2step} performs forward model choice solely on adjusted
+ \eqn{R^2}{R2} and P-value, for ordination objects created by \code{\link{rda}} or \code{\link{capscale}}.
}
\usage{
-ordistep(object, scope, direction = c("both", "backward", "forward"), Pin = 0.05, Pout = 0.1, pstep = 100, perm.max = 1000, steps = 50, trace = TRUE, ...)
+ordistep(object, scope, direction = c("both", "backward", "forward"),
+ Pin = 0.05, Pout = 0.1, pstep = 100, perm.max = 1000, steps = 50,
+ trace = TRUE, ...)
+ordiR2step(object, scope, direction = c("both", "forward"),
+ Pin = 0.05, pstep = 100, perm.max = 1000,
+ trace = TRUE, ...)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{object}{
- An ordination object inheriting from \code{\link{cca}}.
+ In \code{ordistep}, an ordination object inheriting from \code{\link{cca}} or \code{\link{rda}}. In \code{ordiR2step}, the object must inherit from \code{\link{rda}}, that is, it must have been computed using \code{\link{rda}} or \code{\link{capscale}}.
}
\item{scope}{
Defines the range of models examined in the stepwise search.
This should be either a single formula, or a list containing
components \code{upper} and \code{lower}, both formulae.
- See \code{\link{step}} for details.
+ See \code{\link{step}} for details. In \code{ordiR2step}, this defines the
+ upper scope; it can also be an ordination object from with the model is extracted.
}
\item{direction}{
The mode of stepwise search, can be one of \code{"both"},
\code{"backward"}, or \code{"forward"}, with a default of \code{"both"}.
- If the \code{scope} argument is missing the default for \code{direction}
- is \code{"backward"}.
+ If the \code{scope} argument is missing, the default for \code{direction}
+ is \code{"backward"}.
}
\item{Pin, Pout}{
Limits of permutation \eqn{P}-values for adding (\code{Pin}) a term to
@@ -75,16 +84,37 @@
is selected adaptively with respect to the defined decision limit. It
is often sensible to have \code{Pout} \eqn{>} \code{Pin} in stepwise
models to avoid cyclic adds and drops of single terms.
+
+ Function \code{ordiR2step} builds model so that it maximizes adjusted
+ \eqn{R^2}{R2} (function \code{\link{RsquareAdj}}) at every step, and
+ stopping when the adjusted \eqn{R^2}{R2} starts to decrease, or the
+ adjusted \eqn{R^2}{R2} of the \code{scope} is exceeded, or the
+ selected permutation \eqn{P}-value is exceeded (Blanchet et
+ al. 2008). The \code{direction} has choices \code{"forward"} and
+ \code{"both"}, but it is very excepctional that a term is dropped with
+ the adjusted \eqn{R^2}{R2} criterion. Function uses adjusted
+ \eqn{R^2}{R2} as the criterion, and it cannot be used if the criterion
+ cannot be calculated. Therefore it is unavailable for
+ \code{\link{cca}}.
+
+ Functions \code{ordistep} (based on \eqn{P} values) and \code{ordiR2step}
+ (based on adjusted \eqn{R^2}{R2} and hence on eigenvalues) can select
+ variables in different order.
}
\value{
- Function returns the selected model with one additional
+ Functions return the selected model with one additional
component, \code{anova}, which contains brief information of steps
taken. You can suppress voluminous output during model building by
setting \code{trace = FALSE}, and find the summary of model history
in the \code{anova} item.
}
+\references{
+ Blanchet, F. G., Legendre, P. & Borcard, D. (2008) Forward selection
+ of explanatory variables. \emph{Ecology} 89, 2623--2632.
+}
+
\author{
Jari Oksanen
}
@@ -95,15 +125,44 @@
are \code{\link{add1.cca}} and \code{\link{drop1.cca}}, and the
function is modelled after standard \code{\link{step}} (which also can
be used directly but uses AIC for model choice, see
- \code{\link{extractAIC.cca}}).
+ \code{\link{extractAIC.cca}}). Function \code{ordiR2step} builds upon
+ \code{\link{RsquareAdj}}.
}
\examples{
## See add1.cca for another example
+
+### Dune data
data(dune)
data(dune.env)
-mod1 <- rda(dune ~ ., dune.env)
-ordistep(mod1, perm.max = 200)
+mod0 <- rda(dune ~ 1, dune.env) # Model with intercept only
+mod1 <- rda(dune ~ ., dune.env) # Model with all explanatory variables
+
+## With scope present, the default direction is "both"
+ordistep(mod0, scope = formula(mod1), perm.max = 200)
ordistep(rda(dune ~ 1, dune.env), scope = formula(mod1), perm.max = 200)
+
+## Example without scope. Default direction is "backward"
+ordistep(mod1, perm.max = 200)
+
+## Example of ordistep, forward
+ordistep(mod0, scope = formula(mod1), direction="forward", perm.max = 200)
+
+### Mite data
+data(mite)
+data(mite.env)
+mite.hel = decostand(mite, "hel")
+mod0 <- rda(mite.hel ~ 1, mite.env) # Model with intercept only
+mod1 <- rda(mite.hel ~ ., mite.env) # Model with all explanatory variables
+
+## Example of ordiR2step with default direction = "both"
+## (This never goes "backward" but evaluates included terms.)
+step.res <- ordiR2step(mod0, mod1, perm.max = 200)
+step.res$anova # Summary table
+
+## Example of ordiR2step with direction = "forward"
+step.res <- ordiR2step(mod0, scope = formula(mod1), direction="forward")
+step.res <- ordiR2step(mod0, scope = formula(mod1), direction="forward", trace=0)
+step.res$anova # Summary table
}
\keyword{ multivariate }
More information about the Vegan-commits
mailing list