[Vegan-commits] r1623 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 4 12:19:15 CEST 2011
Author: jarioksa
Date: 2011-06-04 12:19:14 +0200 (Sat, 04 Jun 2011)
New Revision: 1623
Modified:
branches/1.17/R/capscale.R
branches/1.17/R/cca.default.R
branches/1.17/R/cca.formula.R
branches/1.17/R/print.cca.R
branches/1.17/R/rda.default.R
branches/1.17/R/rda.formula.R
branches/1.17/inst/ChangeLog
branches/1.17/man/cca.object.Rd
Log:
merged capscale upgrade for missing 'comm' & cca/rda/capscale handling of zeroed components
Modified: branches/1.17/R/capscale.R
===================================================================
--- branches/1.17/R/capscale.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/capscale.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -13,8 +13,10 @@
data <- ordiGetData(match.call(), environment(formula))
}
formula <- formula(terms(formula, data = data))
- X <- formula[[2]]
- X <- eval(X, environment(formula))
+ ## The following line was eval'ed in environment(formula), but
+ ## that made update() fail. Rethink the line if capscale() fails
+ ## mysteriously at this point.
+ X <- eval(formula[[2]]) #, environment(formula))
if (!inherits(X, "dist")) {
comm <- X
dfun <- match.fun(dfun)
@@ -93,19 +95,19 @@
colnames(sol$CCA$wa) <- colnames(sol$CCA$v) <-
paste("CAP", 1:ncol(sol$CCA$u), sep = "")
}
- if (!is.null(sol$CA)) {
+ if (!is.null(sol$CA) && sol$CA$rank > 0) {
colnames(sol$CA$u) <- names(sol$CA$eig) <- colnames(sol$CA$v) <-
paste("MDS", 1:ncol(sol$CA$u), sep = "")
- ## update for negative eigenvalues
- poseig <- length(sol$CA$eig)
- if (any(X$eig < 0)) {
- negax <- X$eig[X$eig < 0]
- sol$CA$imaginary.chi <- sum(negax)
- sol$tot.chi <- sol$tot.chi + sol$CA$imaginary.chi
- sol$CA$imaginary.rank <- length(negax)
- sol$CA$imaginary.u.eig <- X$negaxes
- }
}
+ ## update for negative eigenvalues
+ poseig <- length(sol$CA$eig)
+ if (any(X$eig < 0)) {
+ negax <- X$eig[X$eig < 0]
+ sol$CA$imaginary.chi <- sum(negax)
+ sol$tot.chi <- sol$tot.chi + sol$CA$imaginary.chi
+ sol$CA$imaginary.rank <- length(negax)
+ sol$CA$imaginary.u.eig <- X$negaxes
+ }
if (!is.null(comm)) {
comm <- scale(comm, center = TRUE, scale = FALSE)
sol$colsum <- sd(comm)
@@ -123,11 +125,18 @@
"/")
comm <- qr.resid(sol$CCA$QR, comm)
}
- if (!is.null(sol$CA)) {
+ if (!is.null(sol$CA) && sol$CA$rank > 0) {
sol$CA$v.eig <- t(comm) %*% sol$CA$u/sqrt(k)
sol$CA$v <- sweep(sol$CA$v.eig, 2, sqrt(sol$CA$eig[1:poseig]),
"/")
}
+ } else {
+ ## input data were dissimilarities, and no 'comm' defined:
+ ## species scores make no sense and are made NA
+ sol$CA$v.eig[] <- sol$CA$v[] <- NA
+ if (!is.null(sol$CCA))
+ sol$CCA$v.eig[] <- sol$CCA$v[] <- NA
+ sol$colsum <- NA
}
if (!is.null(sol$CCA))
sol$CCA$centroids <- centroids.cca(sol$CCA$wa, d$modelframe)
Modified: branches/1.17/R/cca.default.R
===================================================================
--- branches/1.17/R/cca.default.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/cca.default.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -39,11 +39,12 @@
Fit = Z, envcentre = attr(Z.r, "centre"))
Xbar <- qr.resid(Q, Xbar)
}
+ if (tmp < ZERO)
+ pCCA$tot.chi <- 0
}
else Z.r <- NULL
if (!missing(Y) && !is.null(Y)) {
Y <- as.matrix(Y)
- rawmat <- Y
Y.r <- weight.centre(Y, rowsum)
Q <- qr(cbind(Z.r, Y.r), tol = ZERO)
if (is.null(pCCA))
@@ -86,42 +87,59 @@
CCA$QR <- Q
CCA$envcentre <- attr(Y.r, "centre")
CCA$Xbar <- Xbar
- Xbar <- qr.resid(Q, Xbar)
- if (exists("exclude.spec")) {
- attr(CCA$v, "na.action") <- exclude.spec
- attr(CCA$v.eig, "na.action") <- exclude.spec
- }
+ } else { # zero rank
+ CCA <- list(eig = 0, rank = rank, qrank = qrank, tot.chi = 0,
+ QR = Q, Xbar = Xbar)
+ u <- matrix(0, nrow=nrow(sol$u), ncol=0)
+ v <- matrix(0, nrow=nrow(sol$v), ncol=0)
+ CCA$u <- CCA$u.eig <- CCA$wa <- CCA$wa.eig <- u
+ CCA$v <- CCA$v.eig <- v
+ CCA$biplot <- matrix(0, 0, 0)
+ CCA$alias <- colnames(Y.r)
}
+ Xbar <- qr.resid(Q, Xbar)
+ if (exists("exclude.spec")) {
+ attr(CCA$v, "na.action") <- exclude.spec
+ attr(CCA$v.eig, "na.action") <- exclude.spec
+ }
+
}
Q <- qr(Xbar)
- if (Q$rank > 0) {
- sol <- svd(Xbar)
- ax.names <- paste("CA", 1:length(sol$d), sep = "")
- colnames(sol$u) <- ax.names
- colnames(sol$v) <- ax.names
- names(sol$d) <- ax.names
- rownames(sol$u) <- rownames(X)
- rownames(sol$v) <- colnames(X)
- rank <- min(Q$rank, sum(sol$d > ZERO))
- if (rank) {
- CA <- list(eig = sol$d[1:rank]^2)
- CA$u <- sweep(as.matrix(sol$u[, 1:rank, drop = FALSE]),
- 1, 1/sqrt(rowsum), "*")
- CA$v <- sweep(as.matrix(sol$v[, 1:rank, drop = FALSE]),
- 1, 1/sqrt(colsum), "*")
- CA$u.eig <- sweep(CA$u, 2, sol$d[1:rank], "*")
- CA$v.eig <- sweep(CA$v, 2, sol$d[1:rank], "*")
- CA$rank <- rank
- CA$tot.chi <- sum(CA$eig)
- CA$Xbar <- Xbar
- if (exists("exclude.spec")) {
- attr(CA$v, "na.action") <- exclude.spec
- attr(CA$v.eig, "na.action") <- exclude.spec
- }
- }
+ sol <- svd(Xbar)
+ ax.names <- paste("CA", 1:length(sol$d), sep = "")
+ colnames(sol$u) <- ax.names
+ colnames(sol$v) <- ax.names
+ names(sol$d) <- ax.names
+ rownames(sol$u) <- rownames(X)
+ rownames(sol$v) <- colnames(X)
+ rank <- min(Q$rank, sum(sol$d > ZERO))
+ if (rank) {
+ CA <- list(eig = sol$d[1:rank]^2)
+ CA$u <- sweep(as.matrix(sol$u[, 1:rank, drop = FALSE]),
+ 1, 1/sqrt(rowsum), "*")
+ CA$v <- sweep(as.matrix(sol$v[, 1:rank, drop = FALSE]),
+ 1, 1/sqrt(colsum), "*")
+ CA$u.eig <- sweep(CA$u, 2, sol$d[1:rank], "*")
+ CA$v.eig <- sweep(CA$v, 2, sol$d[1:rank], "*")
+ CA$rank <- rank
+ CA$tot.chi <- sum(CA$eig)
+ CA$Xbar <- Xbar
+
+ } else { # zero rank: no residual component
+ CA <- list(eig = 0, rank = rank, tot.chi = 0,
+ Xbar = Xbar)
+ CA$u <- CA$u.eig <- matrix(0, nrow(sol$u), 0)
+ CA$v <- CA$v.eig <- matrix(0, nrow(sol$v), 0)
}
+ if (exists("exclude.spec")) {
+ attr(CA$v, "na.action") <- exclude.spec
+ attr(CA$v.eig, "na.action") <- exclude.spec
+ }
call <- match.call()
call[[1]] <- as.name("cca")
+ ## computed pCCA$rank was needed before, but zero it here
+ if (!is.null(pCCA) && pCCA$tot.chi == 0)
+ pCCA$rank <- 0
sol <- list(call = call, grand.total = gran.tot, rowsum = rowsum,
colsum = colsum, tot.chi = tot.chi, pCCA = pCCA, CCA = CCA,
CA = CA)
Modified: branches/1.17/R/cca.formula.R
===================================================================
--- branches/1.17/R/cca.formula.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/cca.formula.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -9,7 +9,7 @@
d <- ordiParseFormula(formula, data, na.action = na.action,
subset = substitute(subset))
sol <- cca.default(d$X, d$Y, d$Z)
- if (!is.null(sol$CCA))
+ if (!is.null(sol$CCA) && sol$CCA$rank > 0)
sol$CCA$centroids <- centroids.cca(sol$CCA$wa, d$modelframe,
sol$rowsum)
if (!is.null(sol$CCA$alias))
Modified: branches/1.17/R/print.cca.R
===================================================================
--- branches/1.17/R/print.cca.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/print.cca.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -24,6 +24,9 @@
rownames(tbl) <- rn[c(TRUE, !is.null(x$CA$imaginary.chi), !is.null(x$pCCA),
!is.null(x$CCA), !is.null(x$CA),
!is.null(x$CA$imaginary.chi))]
+ ## Remove "Proportion" if only one component
+ if (is.null(x$CCA) && is.null(x$pCCA))
+ tbl <- tbl[,-2]
printCoefmat(tbl, digits = digits, na.print = "")
cat("Inertia is", x$inertia, "\n")
if (!is.null(x$CCA$alias))
@@ -37,11 +40,11 @@
cat(length(sp.na), "species",
ifelse(length(sp.na)==1, "(variable)", "(variables)"),
"deleted due to missingness\n")
- if (!is.null(x$CCA)) {
+ if (!is.null(x$CCA) && x$CCA$rank > 0) {
cat("\nEigenvalues for constrained axes:\n")
print(x$CCA$eig, digits = digits, ...)
}
- if (!is.null(x$CA)) {
+ if (!is.null(x$CA) && x$CA$rank > 0) {
ax.lim <- 8
ax.trig <- 16
cat("\nEigenvalues for unconstrained axes:\n")
Modified: branches/1.17/R/rda.default.R
===================================================================
--- branches/1.17/R/rda.default.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/rda.default.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -23,11 +23,12 @@
Fit = Z, envcentre = attr(Z.r, "scaled:center"))
Xbar <- qr.resid(Q, Xbar)
}
+ if (tmp < ZERO)
+ pCCA$tot.chi <- 0
}
else Z.r <- NULL
if (!missing(Y) && !is.null(Y)) {
Y <- as.matrix(Y)
- rawmat <- Y
Y.r <- scale(Y, center = TRUE, scale = FALSE)
Q <- qr(cbind(Z.r, Y.r), tol = ZERO)
if (is.null(pCCA))
@@ -72,31 +73,43 @@
CCA$envcentre <- attr(Y.r, "scaled:center")
CCA$Xbar <- Xbar
Xbar <- qr.resid(Q, Xbar)
+ } else {
+ CCA <- list(eig = 0, rank = rank, qrank = qrank, tot.chi = 0,
+ QR = Q, Xbar = Xbar)
+ u <- matrix(0, nrow=nrow(sol$u), ncol=0)
+ v <- matrix(0, nrow=nrow(sol$v), ncol=0)
+ CCA$u <- CCA$u.eig <- CCA$wa <- CCA$wa.eig <- u
+ CCA$v <- CCA$v.eig <- v
+ CCA$biplot <- matrix(0, 0, 0)
+ CCA$alias <- colnames(Y.r)
}
}
Q <- qr(Xbar)
- if (Q$rank > 0) {
- sol <- svd(Xbar)
- sol$d <- sol$d/sqrt(NR)
- ax.names <- paste("PC", 1:length(sol$d), sep = "")
- colnames(sol$u) <- ax.names
- colnames(sol$v) <- ax.names
- names(sol$d) <- ax.names
- rownames(sol$u) <- rownames(X)
- rownames(sol$v) <- colnames(X)
- rank <- min(Q$rank, sum(sol$d > ZERO))
- if (rank) {
- CA <- list(eig = (sol$d[1:rank]^2))
- CA$u <- as.matrix(sol$u)[, 1:rank, drop = FALSE]
- CA$v <- as.matrix(sol$v)[, 1:rank, drop = FALSE]
- CA$u.eig <- sweep(as.matrix(CA$u), 2, sol$d[1:rank],
- "*")
- CA$v.eig <- sweep(as.matrix(CA$v), 2, sol$d[1:rank],
- "*")
- CA$rank <- rank
- CA$tot.chi <- sum(CA$eig)
- CA$Xbar <- Xbar
- }
+ sol <- svd(Xbar)
+ sol$d <- sol$d/sqrt(NR)
+ ax.names <- paste("PC", 1:length(sol$d), sep = "")
+ colnames(sol$u) <- ax.names
+ colnames(sol$v) <- ax.names
+ names(sol$d) <- ax.names
+ rownames(sol$u) <- rownames(X)
+ rownames(sol$v) <- colnames(X)
+ rank <- min(Q$rank, sum(sol$d > ZERO))
+ if (rank) {
+ CA <- list(eig = (sol$d[1:rank]^2))
+ CA$u <- as.matrix(sol$u)[, 1:rank, drop = FALSE]
+ CA$v <- as.matrix(sol$v)[, 1:rank, drop = FALSE]
+ CA$u.eig <- sweep(as.matrix(CA$u), 2, sol$d[1:rank],
+ "*")
+ CA$v.eig <- sweep(as.matrix(CA$v), 2, sol$d[1:rank],
+ "*")
+ CA$rank <- rank
+ CA$tot.chi <- sum(CA$eig)
+ CA$Xbar <- Xbar
+ } else { # zero rank: no residual component
+ CA <- list(eig = 0, rank = rank, tot.chi = 0,
+ Xbar = Xbar)
+ CA$u <- CA$u.eig <- matrix(0, nrow(sol$u), 0)
+ CA$v <- CA$v.eig <- matrix(0, nrow(sol$v), 0)
}
call <- match.call()
call[[1]] <- as.name("rda")
Modified: branches/1.17/R/rda.formula.R
===================================================================
--- branches/1.17/R/rda.formula.R 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/R/rda.formula.R 2011-06-04 10:19:14 UTC (rev 1623)
@@ -1,6 +1,6 @@
-"rda.formula" <-
-function (formula, data, scale = FALSE, na.action = na.fail,
- subset = NULL, ...)
+`rda.formula` <-
+ function (formula, data, scale = FALSE, na.action = na.fail,
+ subset = NULL, ...)
{
if (missing(data)) {
data <- parent.frame()
@@ -10,7 +10,7 @@
d <- ordiParseFormula(formula, data, na.action = na.action,
subset = substitute(subset))
sol <- rda.default(d$X, d$Y, d$Z, scale)
- if (!is.null(sol$CCA))
+ if (!is.null(sol$CCA) && sol$CCA$rank > 0)
sol$CCA$centroids <- centroids.cca(sol$CCA$wa, d$modelframe)
if (!is.null(sol$CCA$alias))
sol$CCA$centroids <- unique(sol$CCA$centroids)
Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/inst/ChangeLog 2011-06-04 10:19:14 UTC (rev 1623)
@@ -4,12 +4,23 @@
Version 1.17-11 (opened April 29, 2011)
+ * merged r1622: document zeroed components in cca.object.Rd.
+
+ * merged r1620: print.cca does not show "Proportions" in
+ unconstrained models with only one component.
+
+ * merged r1617 thru 1619: zero components in rda() and capscale()
+ similarly as r1517 for cca().
+
* merged r1615: inconsistent notatation in MOStest.Rd.
* merged r1614: ordiR2step more informative.
* merged r1609: specnumber gained 'groups'.
+ * merged r1608, 10, 11, 13: capscale() sets species scores to NA
+ if input were dissimilarities and no 'comm' was given.
+
* merged r1603: cca/rda support functions failed with prc() and
now refuse its results.
@@ -23,6 +34,8 @@
* merged r1588: speed-up of examples.
+ * merged r1517: zero components in cca().
+
* nobs: copied nobs() methods for vegan from r1596.
Version 1.17-10 (released April 27, 2011)
Modified: branches/1.17/man/cca.object.Rd
===================================================================
--- branches/1.17/man/cca.object.Rd 2011-06-04 09:56:02 UTC (rev 1622)
+++ branches/1.17/man/cca.object.Rd 2011-06-04 10:19:14 UTC (rev 1623)
@@ -46,13 +46,22 @@
values). This is a vector of indices of missing value rows in the
original data and a class of the action, usually either
\code{"omit"} or \code{"exclude"}.}
-\item{pCCA, CCA, CA}{Actual ordination results for conditioned
+
+ \item{pCCA, CCA, CA}{Actual ordination results for conditioned
(partial), constrained and unconstrained components of the
- model. Any of these can be \code{NULL} if there is no corresponding
- component.
- Items \code{pCCA}, \code{CCA} and \code{CA} have similar
- structure, and contain following items:
- \describe{
+ model. If constraints or conditions are not given, the
+ corresponding components \code{CCA} and \code{pCCA} are
+ \code{NULL}. If they are specified but have zero rank and zero
+ eigenvalue (e.g., due to aliasing), they have a standard structure
+ like described below, but the result scores have zero columns, but
+ the correct number of rows. The residual component is never
+ \code{NULL}, and if there is no residual variation (like in
+ overdefined model), its scores have zero columns. The standard
+ \code{print} command does not show \code{NULL} components, but it
+ prints zeros for zeroed components. Items \code{pCCA}, \code{CCA}
+ and \code{CA} contain following items:
+
+\describe{
\item{\code{alias}}{The names of the aliased constraints or conditions.
Function \code{\link{alias.cca}} does not access this item
directly, but it finds the aliased variables and their defining
More information about the Vegan-commits
mailing list