[Vegan-commits] r853 - in branches/1.15: R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Jun 2 18:57:07 CEST 2009


Author: jarioksa
Date: 2009-06-02 18:57:07 +0200 (Tue, 02 Jun 2009)
New Revision: 853

Added:
   branches/1.15/R/eigenvals.R
   branches/1.15/R/simulate.rda.R
   branches/1.15/man/eigenvals.Rd
   branches/1.15/man/simulate.rda.Rd
Modified:
   branches/1.15/R/print.summary.cca.R
   branches/1.15/R/summary.cca.R
   branches/1.15/inst/ChangeLog
Log:
copied simulate.rda and eigenvals and merged r845 (summary.cca uses eigenvals) to branches/1.15

Copied: branches/1.15/R/eigenvals.R (from rev 852, pkg/vegan/R/eigenvals.R)
===================================================================
--- branches/1.15/R/eigenvals.R	                        (rev 0)
+++ branches/1.15/R/eigenvals.R	2009-06-02 16:57:07 UTC (rev 853)
@@ -0,0 +1,85 @@
+# Extract eigenvalues from an object that has them
+
+`eigenvals` <-
+    function(x, ...)
+{
+    UseMethod("eigenvals")
+}
+
+`eigenvals.default`<-
+    function(x, ...)
+{
+    ## svd and eigen return unspecified 'list', see if this could be
+    ## either of them
+    out <- NA
+    if (is.list(x)) {
+        ## eigen
+        if (length(x) == 2 && all(names(x) %in% c("values", "vectors")))
+            out <- x$values
+        ## svd: return squares of singular values
+        if (length(x) == 3 && all(names(x) %in% c("d", "u", "v")))
+            out <- x$d^2
+    }
+    class(out) <- "eigenvals"
+    out
+}
+
+## squares of sdev 
+`eigenvals.prcomp` <-
+    function(x, ...)
+{
+    out <- x$sdev^2
+    names(out) <- colnames(x$rotation)
+    class(out) <- "eigenvals"
+    out
+}
+
+## squares of sdev
+`eigenvals.princomp` <-
+    function(x, ...)
+{
+    out <- x$sdev^2
+    class(out) <- "eigenvals"
+    out
+}
+
+## concatenate constrained and unconstrained eigenvalues in cca, rda
+## and capscale (vegan) -- ignore pCCA component
+`eigenvals.cca` <- function(x, constrained = FALSE, ...)
+{
+   if (constrained)
+       out <- x$CCA$eig
+   else
+       out <- c(x$CCA$eig, x$CA$eig)
+   if (!is.null(out))
+       class(out) <- c("eigenvals")
+   out
+}
+
+## wcmdscale (in vegan)
+`eigenvals.wcmdscale` <-
+    function(x, ...)
+{
+    out <- x$eig
+    class(out) <-"eigenvals"
+    out
+}
+
+`print.eigenvals` <-
+    function(x, ...)
+{
+    print(unclass(x), ...)
+    invisible(x)
+}
+
+`summary.eigenvals` <-
+    function(object, ...)
+{
+    vars <- object/sum(object)
+    importance <- rbind(`Eigenvalue` = object,
+                        `Proportion Explained` = round(vars, 5),
+                        `Cumulative Proportion`= round(cumsum(vars), 5))
+    out <- list(importance = importance)
+    class(out) <- c("summary.eigenvals", "summary.prcomp")
+    out
+}

Modified: branches/1.15/R/print.summary.cca.R
===================================================================
--- branches/1.15/R/print.summary.cca.R	2009-06-02 16:35:56 UTC (rev 852)
+++ branches/1.15/R/print.summary.cca.R	2009-06-02 16:57:07 UTC (rev 853)
@@ -22,12 +22,10 @@
         cat("after removing the contribution of conditiniong variables\n")
     }
     cat("\n")
-    out <- rbind("Eig.value" = c(x$ev.con, x$ev.uncon),
-                 "Accounted" = c(x$ev.con.account, x$ev.uncon.account))
-    print(out, digits = digits, ...)
-    if (!is.null(x$cca.acc)) {
+    print(x$cont, ...)
+    if (!is.null(x$concont)) {
         cat("\nAccumulated constrained eigenvalues\n")
-        print(x$cca.acc, digits = digits, ...)
+        print(x$concont, ...)
     }
     cat("\nScaling", x$scaling, "for species and site scores\n")
     if (abs(x$scaling) == 2) {

Copied: branches/1.15/R/simulate.rda.R (from rev 852, pkg/vegan/R/simulate.rda.R)
===================================================================
--- branches/1.15/R/simulate.rda.R	                        (rev 0)
+++ branches/1.15/R/simulate.rda.R	2009-06-02 16:57:07 UTC (rev 853)
@@ -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: branches/1.15/R/summary.cca.R
===================================================================
--- branches/1.15/R/summary.cca.R	2009-06-02 16:35:56 UTC (rev 852)
+++ branches/1.15/R/summary.cca.R	2009-06-02 16:57:07 UTC (rev 853)
@@ -23,17 +23,8 @@
     summ$partial.chi <- object$pCCA$tot.chi
     summ$constr.chi <- object$CCA$tot.chi
     summ$unconst.chi <- object$CA$tot.chi
-    summ$ev.con <- object$CCA$eig
-    summ$ev.uncon <- object$CA$eig
-    ev.account <- summ$tot.chi
-    if (!is.null(object$pCCA)) 
-        ev.account <- ev.account - summ$partial.chi
-    if (!is.null(object$CCA))
-        summ$ev.con.account <- cumsum(summ$ev.con)/ev.account
-    summ$ev.uncon.account <-
-        (max(summ$constr.chi, 0) + cumsum(summ$ev.uncon))/ev.account
-    if (!is.null(object$CCA))
-        summ$cca.acc <- cumsum(summ$ev.con)/summ$constr.chi
+    summ$cont <- summary(eigenvals(object))
+    summ$concont <- summary(eigenvals(object, constrained = TRUE))
     summ$ev.head <- c(summ$ev.con, summ$ev.uncon)[1:axes]
     summ$scaling <- scaling
     summ$digits <- digits

Modified: branches/1.15/inst/ChangeLog
===================================================================
--- branches/1.15/inst/ChangeLog	2009-06-02 16:35:56 UTC (rev 852)
+++ branches/1.15/inst/ChangeLog	2009-06-02 16:57:07 UTC (rev 853)
@@ -5,6 +5,12 @@
 
 Version 1.15-3 (opened 28 May, 2009)
 
+	* merged r845: summary.cca uses eigenvals.
+	
+	* copied eigenvals at r845.
+	
+	* copied simulate.rda at r814.
+
 	* merged r818, 820: tsallis got argument hill.
 
 	* merged r796, 801: ordisurf can do linear or quadratic surfaces

Copied: branches/1.15/man/eigenvals.Rd (from rev 852, pkg/vegan/man/eigenvals.Rd)
===================================================================
--- branches/1.15/man/eigenvals.Rd	                        (rev 0)
+++ branches/1.15/man/eigenvals.Rd	2009-06-02 16:57:07 UTC (rev 853)
@@ -0,0 +1,84 @@
+\name{eigenvals}
+\Rdversion{1.1}
+\alias{eigenvals}
+\alias{eigenvals.default}
+\alias{eigenvals.prcomp}
+\alias{eigenvals.princomp}
+\alias{eigenvals.cca}
+\alias{eigenvals.wcmdscale}
+\alias{print.eigenvals}
+\alias{summary.eigenvals}
+
+\title{
+  Extract Eigenvalues from an Ordination Object
+}
+\description{
+  Function extracts eigenvalues from an object that has them. Many
+  multivariate methods return such objects. 
+}
+\usage{
+eigenvals(x, ...)
+\method{eigenvals}{cca}(x, constrained = FALSE, ...)
+\method{summary}{eigenvals}(object, ...)
+}
+
+\arguments{
+  \item{x}{
+    An object from which to extract eigenvalues.
+}
+  \item{object}{
+    An \code{eigenvals} result object.
+}
+   \item{constrained}{
+    Return only constrained eigenvalues.
+}
+   \item{\dots}{
+    Other arguments to the functions (usually ignored)
+}
+}
+
+\details{
+  This is a generic function that has methods for \code{\link{cca}},
+  \code{\link{wcmdscale}}, \code{\link{prcomp}} and
+  \code{\link{princomp}} result objects. The default method also
+  extracts eigenvalues if the result looks like being from
+  \code{\link{eigen}} or \code{\link{svd}}.  Functions
+  \code{\link{prcomp}} and \code{\link{princomp}} contain square roots
+  of eigenvalues that all called standard deviations, but
+  \code{eigenvals} function returns their squares.  Function
+  \code{\link{svd}} contains singular values, but function
+  \code{eigenvals} returns their squares. For constrained ordination
+  methods \code{\link{cca}}, \code{\link{rda}} and
+  \code{\link{capscale}} the function returns the both constrained and
+  unconstrained eigenvalues concatenated in one vector, but the partial
+  component will be ignored. However, with argument 
+  \code{constrained = TRUE} only constrained eigenvalues are returned. 
+
+  The \code{summary} of \code{eigenvals} result returns eigenvalues,
+  proportion explained and cumulative proportion explained.
+}
+
+\value{
+  An object of class \code{"eigenvals"} which is a vector of eigenvalues.
+}
+
+\author{
+ Jari Oksanen.
+}
+
+\seealso{
+ \code{\link{eigen}}, \code{\link{svd}}, \code{\link{prcomp}},
+ \code{\link{princomp}}, \code{\link{cca}}, \code{\link{rda}},
+ \code{\link{capscale}}, \code{\link{wcmdscale}},
+ \code{\link{cca.object}}. 
+}
+\examples{
+data(varespec)
+data(varechem)
+mod <- cca(varespec ~ Al + P + K, varechem)
+ev <- eigenvals(mod)
+ev
+summary(ev)
+}
+\keyword{ multivariate }
+

Copied: branches/1.15/man/simulate.rda.Rd (from rev 852, pkg/vegan/man/simulate.rda.Rd)
===================================================================
--- branches/1.15/man/simulate.rda.Rd	                        (rev 0)
+++ branches/1.15/man/simulate.rda.Rd	2009-06-02 16:57:07 UTC (rev 853)
@@ -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 }



More information about the Vegan-commits mailing list