[Vegan-commits] r402 - in pkg: R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Jun 8 18:06:53 CEST 2008


Author: jarioksa
Date: 2008-06-08 18:06:53 +0200 (Sun, 08 Jun 2008)
New Revision: 402

Modified:
   pkg/R/alias.cca.R
   pkg/R/calibrate.cca.R
   pkg/R/cca.default.R
   pkg/R/extractAIC.cca.R
   pkg/R/permutest.cca.R
   pkg/R/predict.cca.R
   pkg/R/predict.rda.R
   pkg/R/print.cca.R
   pkg/R/rda.default.R
   pkg/inst/ChangeLog
   pkg/man/cca.object.Rd
Log:
fixes to rank issues in cca/rda

Modified: pkg/R/alias.cca.R
===================================================================
--- pkg/R/alias.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/alias.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -14,7 +14,7 @@
     R <- R[1:min(dim(R)), , drop = FALSE]
     R[lower.tri(R)] <- 0
     d <- dim(R)
-    rank <- object$CCA$rank
+    rank <- object$CCA$QR$rank
     p <- d[2]
     value$Complete <- if (is.null(p) || rank == p) 
         NULL

Modified: pkg/R/calibrate.cca.R
===================================================================
--- pkg/R/calibrate.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/calibrate.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -5,6 +5,8 @@
         stop("does not work with conditioned (partial) models")
     if (is.null(object$CCA))
         stop("needs constrained model")
+    if (object$CCA$rank < object$CCA$qrank)
+        stop("rank of constraints is higher than rank of dependent data")
     if (rank != "full")
         rank <- min(rank, object$CCA$rank)
     else
@@ -13,11 +15,12 @@
         wa <- object$CCA$wa        
     else
         wa <- predict(object, type="wa", newdata=newdata)
-    b <- (coef(object))[object$CCA$QR$pivot[1:object$CCA$rank], , drop=FALSE]
+    qrank <- object$CCA$qrank
+    b <- (coef(object))[object$CCA$QR$pivot[1:qrank], , drop=FALSE]
     b <- solve(b)
-    pred <- wa[ , 1:rank, drop=FALSE]  %*% b[1:rank, , drop =FALSE]
+    pred <- wa[ , 1:rank, drop=FALSE]  %*% b[1:qrank, , drop =FALSE]
     envcen <- object$CCA$envcentre[object$CCA$QR$pivot]
-    envcen <- envcen[1:object$CCA$rank]
+    envcen <- envcen[1:object$CCA$qrank]
     pred <- sweep(pred, 2, envcen, "+")
     pred
 }

Modified: pkg/R/cca.default.R
===================================================================
--- pkg/R/cca.default.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/cca.default.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -84,6 +84,7 @@
             CCA$biplot <- cor(Y.r[, oo, drop = FALSE], sol$u[, 
                                         1:rank, drop = FALSE])
             CCA$rank <- rank
+            CCA$qrank <- qrank
             CCA$tot.chi <- sum(CCA$eig)
             CCA$QR <- Q
             CCA$envcentre <- attr(Y.r, "centre")

Modified: pkg/R/extractAIC.cca.R
===================================================================
--- pkg/R/extractAIC.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/extractAIC.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -3,7 +3,7 @@
 {
    n <- nrow(fit$CA$Xbar)
    edf <- 1
-   if (!is.null(fit$CCA$rank)) edf <- edf + fit$CCA$rank
+   if (!is.null(fit$CCA$rank)) edf <- edf + fit$CCA$qrank
    if (!is.null(fit$pCCA$rank)) edf <- edf + fit$pCCA$rank
    #edf <- n - fit$CA$rank
    RSS <- deviance(fit)

Modified: pkg/R/permutest.cca.R
===================================================================
--- pkg/R/permutest.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/permutest.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -18,7 +18,7 @@
     else {
         Chi.z <- x$CCA$tot.chi
         names(Chi.z) <- "Model"
-        q <- x$CCA$rank
+        q <- x$CCA$qrank
     }
     Chi.xz <- x$CA$tot.chi
     names(Chi.xz) <- "Residual"

Modified: pkg/R/predict.cca.R
===================================================================
--- pkg/R/predict.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/predict.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -39,7 +39,7 @@
             }
             E <- sweep(E, 2, c(object$pCCA$envcentre, object$CCA$envcentre), 
                        "-")
-            p1 <- object[[model]]$QR$pivot[1:object[[model]]$rank]
+            p1 <- object[[model]]$QR$pivot[1:object[[model]]$qrank]
             u <- E[, p1, drop = FALSE] %*% coef(object)[p1, , 
                          drop = FALSE]
             u <- u[, 1:take, drop = FALSE]

Modified: pkg/R/predict.rda.R
===================================================================
--- pkg/R/predict.rda.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/predict.rda.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -44,7 +44,7 @@
                                       object$terminfo$xlev)
                 E <- cbind(d$Z, d$Y)
             }
-            p1 <- object[[model]]$QR$pivot[1:object[[model]]$rank]
+            p1 <- object[[model]]$QR$pivot[1:object[[model]]$qrank]
             E <- sweep(E, 2, c(object$pCCA$envcentre, object$CCA$envcentre), 
                        "-")
             u <- E[, p1, drop = FALSE] %*% coef(object)[p1, , 

Modified: pkg/R/print.cca.R
===================================================================
--- pkg/R/print.cca.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/print.cca.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -5,12 +5,8 @@
         warning("this is an ade4 object which vegan cannot handle")
         x <- ade2vegancca(x)
     }
-    if (!is.null(x$CCA)) {
-        qrank <- x$CCA$QR$rank
-        if (!is.null(x$pCCA))
-            qrank <- qrank - x$pCCA$rank
-        if (x$CCA$rank < qrank)
-            warning("rank of constraints higher than the rank of dependent data\nvegan may not handle this -- wait for a bug fix release", call. = FALSE)
+    if (!is.null(x$CCA) && x$CCA$rank < x$CCA$qrank) {
+        warning("rank of constraints higher than the rank of dependent data\nvegan may not handle this -- wait for a bug fix release", call. = FALSE)
     }
     cat("\nCall:\n")
     cat(deparse(x$call), "\n\n")

Modified: pkg/R/rda.default.R
===================================================================
--- pkg/R/rda.default.R	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/R/rda.default.R	2008-06-08 16:06:53 UTC (rev 402)
@@ -68,6 +68,7 @@
             CCA$biplot <- cor(Y.r[, oo, drop = FALSE], sol$u[, 
                                         1:rank, drop = FALSE])
             CCA$rank <- rank
+            CCA$qrank <- qrank
             CCA$tot.chi <- sum(CCA$eig)
             CCA$QR <- Q
             CCA$envcentre <- attr(Y.r, "scaled:center")

Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/inst/ChangeLog	2008-06-08 16:06:53 UTC (rev 402)
@@ -14,6 +14,22 @@
 	their rank is higher than the rank of the community data.  Because
 	cca objects are handled by many other functions, this may cause
 	trouble in surprising places and the change is unstabilizing.
+
+	* rda.default, cca.default: return now both 'rank' for the
+	ordination, and 'qrank' or the rank of the constrainst. This means
+	following changes in functions using ranks:
+
+	* alias.cca: qrank
+	* anova.ccabyaxis: OK 
+	* bstick.cca: OK 
+	* calibrate.cca: stop with error if qrank > rank
+	* extractAIC: qrank
+	* mso: OK (HW confirms)
+	* msoplot: OK (HW cofirms)
+	* permutest.cca: qrank, also fixes anova.ccabyterm
+	* predict.cca/rda: qrank with type = "lc"
+	* summary.cca/rda: OK
+	* scores.cca/rda: OK
 	
 Version 1.14-2 (closed June 5, 2008)
 

Modified: pkg/man/cca.object.Rd
===================================================================
--- pkg/man/cca.object.Rd	2008-06-05 10:39:26 UTC (rev 401)
+++ pkg/man/cca.object.Rd	2008-06-08 16:06:53 UTC (rev 402)
@@ -69,7 +69,10 @@
       \code{\link{predict.cca}}, \code{\link{predict.rda}},
       \code{\link{calibrate.cca}}.  For possible uses of this component,
       see \code{\link{qr}}. In \code{pCCA} and \code{CCA}.} 
-    \item{\code{rank}}{The rank of the component.}
+    \item{\code{rank}}{The rank of the ordination component.}
+    \item{\code{qrank}}{The rank of the constraints which is the
+     difference of the ranks of QR decompositions in \code{pCCA} and
+     \code{CCA} components. Only in \code{CCA}.}
     \item{\code{tot.chi}}{Total inertia or the sum of all eigenvalues of the
       component.}
     \item{\code{u}}{(Weighted) orthonormal site scores.  Please note that



More information about the Vegan-commits mailing list