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

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

Author: jarioksa
Date: 2009-11-29 14:57:38 +0100 (Sun, 29 Nov 2009)
New Revision: 1075

Modified:
   pkg/vegan/R/simulate.rda.R
   pkg/vegan/inst/ChangeLog
Log:
implemented simulate.cca

Modified: pkg/vegan/R/simulate.rda.R
===================================================================
--- pkg/vegan/R/simulate.rda.R	2009-11-29 13:06:21 UTC (rev 1074)
+++ pkg/vegan/R/simulate.rda.R	2009-11-29 13:57:38 UTC (rev 1075)
@@ -26,14 +26,54 @@
                nrow = nrow(ftd)))
     else
         ans <- as.data.frame(ftd + object$CA$Xbar[indx,])
-    attr(ans, "seed") <- RNGstate
+    if (is.null(indx))
+        attr(ans, "seed") <- RNGstate
+    else
+        attr(ans, "seed") <- indx
     ans
 }
 
+### simulate. cca was cloned from simulate.rda.  Works with internal
+### Chi-square standardized form, and at the end back-standardizes
+### with row and column totals and matrix grand totals. This does not
+### still guarantee that all marginal totals are positive.
+
 `simulate.cca` <-
-    function(object, nsim = 1, seed = NULL, ...)
+    function(object, nsim = 1, seed = NULL, indx = NULL, ...)
 {
-    .NotYetImplemented()
+    ## 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, type = "working")
+    ## pCCA: add partial Fit to the constrained
+    if (!is.null(object$pCCA))
+        ftd <- ftd + object$pCCA$Fit
+    if (is.null(indx))
+        ans <- as.data.frame(ftd + matrix(rnorm(length(ftd), 
+               sd = outer(rep(1,nrow(ftd)), sd(object$CA$Xbar))), 
+               nrow = nrow(ftd)))
+    else
+        ans <- as.data.frame(ftd + object$CA$Xbar[indx,])
+    ## From internal form to the original form with fixed marginal totals
+    rc <- object$rowsum %o% object$colsum
+    ans <- (ans * sqrt(rc) + rc) * object$grand.total
+    if (is.null(indx))
+        attr(ans, "seed") <- RNGstate
+    else
+        attr(ans, "seed") <- indx
+    ans
 }
 
 `simulate.capscale` <-

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2009-11-29 13:06:21 UTC (rev 1074)
+++ pkg/vegan/inst/ChangeLog	2009-11-29 13:57:38 UTC (rev 1075)
@@ -24,6 +24,12 @@
 	duplicates such as from sample(nrow(data), replace = TRUE) so that
 	bootstrapping is possible unlike in anova.cca/permutest.cca. Works
 	now with partial model: adds partial fit to the constrained fit.
+
+	* simulate.cca: implemented by cloning simulate.rda. Function
+	works with internal Chi-square standardized form and
+	back-transforms the result by marginal totals and matrix grand
+	total at the end. The marginal totals may still be negative. Not
+	explicitly documented, and needs a connoisseur user. 
 	
 Version 1.16-32 (closed November 13, 2009)