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

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

Author: jarioksa
Date: 2009-02-15 13:05:20 +0100 (Sun, 15 Feb 2009)
New Revision: 684

Added:
   pkg/vegan/R/simulate.rda.R
   pkg/vegan/man/simulate.rda.Rd
Modified:
   pkg/vegan/DESCRIPTION
   pkg/vegan/inst/ChangeLog
Log:
added simulate.rda

Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION	2009-02-12 13:06:24 UTC (rev 683)
+++ pkg/vegan/DESCRIPTION	2009-02-15 12:05:20 UTC (rev 684)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 1.16-12
-Date: February 9, 2009
+Version: 1.16-13
+Date: February 15, 2009
 Author: Jari Oksanen, Roeland Kindt, Pierre Legendre, Bob O'Hara, Gavin L. Simpson, 
    Peter Solymos, M. Henry H. Stevens, Helene Wagner  
 Maintainer: Jari Oksanen <jari.oksanen at oulu.fi>

Added: pkg/vegan/R/simulate.rda.R
===================================================================
--- pkg/vegan/R/simulate.rda.R	                        (rev 0)
+++ pkg/vegan/R/simulate.rda.R	2009-02-15 12:05:20 UTC (rev 684)
@@ -0,0 +1,40 @@
+`simulate.rda` <-
+    function(object, nsim = 1, seed = NULL, ...) 
+{
+    ## First check cases that won't work (yet?)
+    if (!is.null(object$pCCA))
+        stop("not yet implemented for partial models")
+    ## Handle RNG: code directly from stats::simulate.lm
+    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
+        runif(1)
+    if (is.null(seed)) 
+        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
+    else {
+        R.seed <- get(".Random.seed", envir = .GlobalEnv)
+        set.seed(seed)
+        RNGstate <- structure(seed, kind = as.list(RNGkind()))
+        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
+    }
+    ## Proper simulation: very similar for simulate.lm, but produces a
+    ## response matrix.
+    if (nsim > 1)
+        .NotYetUsed("nsim")
+    ftd <- fitted(object)
+    ans <- as.data.frame(ftd + matrix(rnorm(length(ftd), 
+        sd = outer(rep(1,nrow(ftd)), sd(object$CA$Xbar))), 
+        nrow = nrow(ftd)))
+    attr(ans, "seed") <- RNGstate
+    ans
+}
+
+`simulate.cca` <-
+    function(object, nsim = 1, seed = NULL, ...)
+{
+    .NotYetImplemented()
+}
+
+`simulate.capscale` <-
+    function(object, nsim = 1, seed = NULL, ...)
+{
+    .NotYetImplemented()
+}

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2009-02-12 13:06:24 UTC (rev 683)
+++ pkg/vegan/inst/ChangeLog	2009-02-15 12:05:20 UTC (rev 684)
@@ -2,8 +2,22 @@
 
 VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
 
-Version 1.6-12 (opened Feb 9, 2009)
+Version 1.6-13 (opened Feb 15, 2009)
 
+	* simulate.rda: new method to simulate responses with random error
+	for an rda() result object. The function uses Normal error, and
+	estimates the sd separately for each species from their rda
+	residuals. Normal error is justified by least squares fitting in
+	rda. This could be made to use permutations of residuals. Not yet
+	implemented for cca or capscale results. For cca we -- at least --
+	need to handle weights, and the capscale result would be rda
+	result, since the method cannot be directly implemented for
+	dissimilarities. The function was triggered by Ben Bolker's email
+	to R-devel mainling list
+	https://stat.ethz.ch/pipermail/r-devel/2009-February/052116.html
+	
+Version 1.6-12 (closed Feb 15, 2009)
+
 	* anova.ccabyaxis: gained keyword 'cutoff' (defaults 1) to break
 	from permutation tests after exceeding the P-value given in the
 	argument. The keyword was introduced because rda and cca are

Added: pkg/vegan/man/simulate.rda.Rd
===================================================================
--- pkg/vegan/man/simulate.rda.Rd	                        (rev 0)
+++ pkg/vegan/man/simulate.rda.Rd	2009-02-15 12:05:20 UTC (rev 684)
@@ -0,0 +1,60 @@
+\name{simulate.rda}
+\alias{simulate.rda}
+\alias{simulate.cca}
+\alias{simulate.capscale}
+\title{ Simulate Responses with Gaussian Error for Redundancy Analysis }
+
+\description{ Function simulates a response data frame so that it adds
+ Gaussian error to the fitted responses of Redundancy Analysis
+ (\code{\link{rda}}). The function is a special case of generic
+ \code{\link{simulate}}, and works similarly as \code{simulate.lm}. 
+}
+
+\usage{
+\method{simulate}{rda}(object, nsim = 1, seed = NULL, ...)
+}
+\arguments{
+  \item{object}{an object representing a fitted \code{\link{rda}} model.}
+  \item{nsim}{number of response vectors to simulate. (Not yet used, and 
+    values above 1 will give an error). }
+  \item{seed}{an object specifying if and how the random number
+    generator should be initialized (\sQuote{seeded}). See 
+    \code{\link{simulate}} for details. }
+  \item{\dots}{additional optional arguments (ignored). }
+}
+
+\details{ The implementation follows \code{"lm"} method of
+  \code{\link{simulate}}, and adds Gaussian (Normal) error to the
+  fitted values (\code{\link{fitted.rda}} using function
+  \code{\link{rnorm}}. The standard deviations are estimated
+  independently for each species (column) from the residuals after
+  fitting the constraints.
+}
+
+\value{ Returns a data frame with similar additional arguments on
+  random number seed as \code{\link{simulate}}.  }
+
+\author{Jari Oksanen}
+
+\note{ The function is not implemented for \code{\link{cca}} or
+  \code{\link{capscale}} objects, but only for \code{\link{rda}}.  
+}
+
+\seealso{ \code{\link{simulate}} for the generic case and for
+  \code{\link{lm}} objects. Function \code{\link{fitted.rda}} returns
+  fitted values without the error component. 
+}
+
+\examples{
+data(dune)
+data(dune.env)
+mod <- rda(dune ~  Moisture + Management, dune.env)
+## One simulation
+update(mod, simulate(mod) ~  .)
+## An impression of confidence regions of site scores
+plot(mod, display="sites")
+for (i in 1:5) lines(procrustes(mod, update(mod, simulate(mod) ~ .)), col="blue")
+}
+\keyword{ models }
+\keyword{ datagen }
+\keyword{ multivariate }