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

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

Author: jarioksa
Date: 2013-09-08 19:38:01 +0200 (Sun, 08 Sep 2013)
New Revision: 2608

Modified:
   pkg/vegan/R/simulate.rda.R
   pkg/vegan/man/simulate.rda.Rd
Log:
untested implementation of correlated species in parametric simulate.rda

Modified: pkg/vegan/R/simulate.rda.R
===================================================================
--- pkg/vegan/R/simulate.rda.R	2013-09-08 15:59:14 UTC (rev 2607)
+++ pkg/vegan/R/simulate.rda.R	2013-09-08 17:38:01 UTC (rev 2608)
@@ -1,6 +1,10 @@
 `simulate.rda` <-
-    function(object, nsim = 1, seed = NULL, indx = NULL, rank = "full", ...) 
+    function(object, nsim = 1, seed = NULL, indx = NULL, rank = "full",
+             iid = TRUE, ...) 
 {
+    ## is.null(indx) && !iid requires MASS
+    if(is.null(indx) && !iid)
+        require(MASS) || stop("simulate options require MASS package")
     ## Handle RNG: code directly from stats::simulate.lm
     if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
         runif(1)
@@ -33,15 +37,27 @@
     ## pRDA: add partial Fit to the constrained
     if (!is.null(object$pCCA))
         ftd <- ftd + object$pCCA$Fit
+    ## if(is.null(indx)), we have parametric Gaussian simulation and
+    ## need to generate sd matrices. The residuals sd is always taken
+    ## from the unconstrained (residual) component $CA$Xbar. If iid,
+    ## we need only species sd's, but if non-independent, we also need
+    ## species covariances.
+    if (iid)
+        dev <- outer(rep(1, nrow(ftd)), apply(object$CA$Xbar, 2, sd))
+    else
+        dev <- cov(object$CA$Xbar)
     ## Generate an array
     ans <- array(0, c(dim(ftd), nsim))
     for (i in seq_len(nsim)) {
-        if (is.null(indx))
-            ans[,,i] <- as.matrix(ftd + matrix(rnorm(length(ftd), 
-               sd = outer(rep(1,nrow(ftd)), apply(object$CA$Xbar, 2, sd))), 
-               nrow = nrow(ftd)))
-        else
+        if (!is.null(indx))
             ans[,,i] <- as.matrix(ftd + object$CA$Xbar[indx[i,],])
+        else if (iid)
+            ans[,,i] <- as.matrix(ftd + matrix(rnorm(length(ftd), sd = dev),
+                                               nrow = nrow(ftd)))
+        else {
+            ans[,,i] <- apply(ftd, 1,
+                              function(x) mvrnorm(1, mu = x, Sigma = dev)) 
+        }
     }
     ## set RNG attributes
     if (is.null(indx))

Modified: pkg/vegan/man/simulate.rda.Rd
===================================================================
--- pkg/vegan/man/simulate.rda.Rd	2013-09-08 15:59:14 UTC (rev 2607)
+++ pkg/vegan/man/simulate.rda.Rd	2013-09-08 17:38:01 UTC (rev 2608)
@@ -12,7 +12,8 @@
  works similarly as \code{simulate.lm}.  }
 
 \usage{
-\method{simulate}{rda}(object, nsim = 1, seed = NULL, indx = NULL, rank = "full", ...)
+\method{simulate}{rda}(object, nsim = 1, seed = NULL, indx = NULL,
+    rank = "full", iid = TRUE, ...)
 }
 \arguments{
 
@@ -40,6 +41,9 @@
 
   \item{rank}{The rank of the constrained component: passed to
     \code{\link{predict.rda}} or \code{\link{predict.cca}}. }
+
+  \item{iid}{Not yet documented (and implementation untested).}
+
   \item{\dots}{additional optional arguments (ignored). }
 }