[Vegan-commits] r1126 - in branches/1.17: . R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Feb 16 07:47:40 CET 2010


Author: jarioksa
Date: 2010-02-16 07:47:38 +0100 (Tue, 16 Feb 2010)
New Revision: 1126

Modified:
   branches/1.17/DESCRIPTION
   branches/1.17/R/CCorA.R
   branches/1.17/R/biplot.CCorA.R
   branches/1.17/R/envfit.default.R
   branches/1.17/R/print.CCorA.R
   branches/1.17/inst/ChangeLog
   branches/1.17/man/
   branches/1.17/man/CCorA.Rd
Log:
merged fixes for CCCorA (r1123) and envfit (r1124) to branches/1.17

Modified: branches/1.17/DESCRIPTION
===================================================================
--- branches/1.17/DESCRIPTION	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/DESCRIPTION	2010-02-16 06:47:38 UTC (rev 1126)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 1.17-0
-Date: January 11, 2010
+Version: 1.17-1
+Date: February 16, 2010
 Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, 
    R. G. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, 
    Helene Wagner  

Modified: branches/1.17/R/CCorA.R
===================================================================
--- branches/1.17/R/CCorA.R	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/R/CCorA.R	2010-02-16 06:47:38 UTC (rev 1126)
@@ -1,19 +1,20 @@
 `CCorA` <-
     function(Y, X, stand.Y = FALSE, stand.X = FALSE, nperm = 0, ...)
 {
-    require(MASS) || stop("requires packages 'MASS'")
+    require(MASS) || stop("Requires package 'MASS'")
     epsilon <- sqrt(.Machine$double.eps)
+    ##
     ## BEGIN: Internal functions
+    ##
     cov.inv <- function(mat, no, epsilon) {
         ## This function returns:
-        ## 1) mat = the matrix of PCA object scores (by SVD);
+        ## 1) mat = matrix F of the principal components (PCA object scores);
         ## 2) S.inv = the inverse of the covariance matrix;
         ## 3) m = the rank of matrix 'mat'
-        ## The inverse of the PCA covariance matrix is simply the
-        ## diagonal matrix of (1/eigenvalues).  If ncol(mat) = 1, the
-        ## inverse of the covariance matrix simply contains
-        ## 1/var(mat).
-        mat <- as.matrix(mat)
+        ## The inverse of the PCA covariance matrix is the diagonal
+        ## matrix of (1/eigenvalues). If ncol(mat) = 1, the
+        ## inverse of the covariance matrix contains 1/var(mat).
+        mat <- as.matrix(mat)    # 'mat' was centred before input to cov.inv
         if(ncol(mat) == 1) {
             S.inv <- as.matrix(1/var(mat))
             m <- 1
@@ -26,7 +27,7 @@
                 m <- mm
             }
             S.inv <- diag(1/S.svd$d[1:m])
-            mat <- mat %*% S.svd$u[,1:m]
+            mat <- mat %*% S.svd$u[,1:m] # S.svd$u = normalized eigenvectors
         }
         list(mat=mat, S.inv=S.inv, m=m)
     }
@@ -56,22 +57,44 @@
         P <- nGE/(nperm+1)
     }
     ## END: internal functions
+    ##
     Y <- as.matrix(Y)
     var.null(Y,1)
     nY <- nrow(Y)
     p <- ncol(Y)
-    Ynoms <- colnames(Y)
+    if(is.null(colnames(Y))) {
+        Ynoms <- paste("VarY", 1:p, sep="")
+        } else {
+        Ynoms <- colnames(Y)
+        }
     X <- as.matrix(X)
     var.null(X,2)
     nX <- nrow(X)
     q <- ncol(X)
-    Xnoms <- colnames(X)
+    if(is.null(colnames(X))) {
+        Xnoms <- paste("VarX", 1:q, sep="")
+        } else {
+        Xnoms <- colnames(X)
+        }
     if(nY != nX) stop("Different numbers of rows in Y and X")
     n <- nY
-    if((p+q) >= (n-1)) stop("Not enough degrees of freedom!")
-    rownoms <- rownames(Y)
+    if((p+q) >= (n-1)) stop("Not enough degrees of freedom: (p+q) >= (n-1)")
+    if(is.null(rownames(X)) & is.null(rownames(Y))) {
+        rownoms <- paste("Obj", 1:n, sep="")
+        } else {
+        if(is.null(rownames(X))) {
+            rownoms <- rownames(Y)
+            } else {
+            rownoms <- rownames(X)
+            }
+        }
     Y.c <- scale(Y, center = TRUE, scale = stand.Y)
     X.c <- scale(X, center = TRUE, scale = stand.X)
+    ## Check for identical matrices
+    if(p == q) {
+    	if(sum(abs(Y-X)) < epsilon^2) stop("Y and X are identical")
+    	if(sum(abs(Y.c-X.c)) < epsilon^2) stop("After centering, Y and X are identical")
+    	}
     ## Replace Y.c and X.c by tables of their PCA object scores, computed by SVD
     temp <- cov.inv(Y.c, 1, epsilon)
     Y <- temp$mat
@@ -79,7 +102,7 @@
     rownames(Y) <- rownoms
     temp <- cov.inv(X.c, 2, epsilon)
     X <- temp$mat
-     qq <- temp$m
+    qq <- temp$m
     rownames(X) <- rownoms
     ## Covariance matrices, etc. from the PCA scores
     S11 <- cov(Y)
@@ -95,7 +118,11 @@
     ## K summarizes the correlation structure between the two sets of variables
     K <- t(S11.chol.inv) %*% S12 %*% S22.chol.inv
     K.svd <- svd(K)
-    EigenValues <- K.svd$d^2
+    Eigenvalues <- K.svd$d^2
+    ##
+    ## Check for circular covariance matrix
+    if((p == q) & (var(K.svd$d) < epsilon))
+    	cat("Warning: [nearly] circular covariance matrix. The solution may be meaningless.",'\n')
     ## K.svd$u %*% diag(K.svd$d) %*% t(K.svd$v)   # To check that K = U D V'
     axenames <- paste("CanAxis",1:length(K.svd$d),sep="")
     U <- K.svd$u
@@ -104,21 +131,19 @@
     B <- S22.chol.inv %*% V
     Cy <- (Y %*% A)/sqrt(n-1)
     Cx <- (X %*% B)/sqrt(n-1)
-    ## Compute the 'Biplot scores of Y variables' a posteriori --
-    ## use 'ginv' for inversion in case there is collinearity
-    ## AA <- coefficients of the regression of Cy on Y.c, times sqrt(n-1)
-    ## AA <- sqrt(n-1) * [Y'Y]-1 Y' Cy
-    YprY <- t(Y.c) %*% Y.c
-    AA <- sqrt(n-1) * ginv(YprY) %*% t(Y.c) %*% Cy
-    ##
-    ## Compute the 'Biplot scores of X variables' a posteriori --
-    XprX <- t(X) %*% X
-    BB <- sqrt(n-1) * ginv(XprX) %*% t(X) %*% Cx
+    ## Compute the 'Biplot scores of Y and X variables' a posteriori --
+    corr.Y.Cy <- cor(Y.c, Cy)  # To plot Y in biplot in space Y
+    corr.Y.Cx <- cor(Y.c, Cx)  # Available for plotting Y in space of X
+    corr.X.Cy <- cor(X.c, Cy)  # Available for plotting X in space of Y
+    corr.X.Cx <- cor(X.c, Cx)  # To plot X in biplot in space X
     ## Add row and column names
-    rownames(AA) <- Ynoms
-    rownames(BB) <- Xnoms
     rownames(Cy) <- rownames(Cx) <- rownoms
-    colnames(AA) <- colnames(BB) <- colnames(Cy) <- colnames(Cx) <- axenames
+	colnames(Cy) <- colnames(Cx) <- axenames    
+    rownames(corr.Y.Cy) <- rownames(corr.Y.Cx) <- Ynoms
+    rownames(corr.X.Cy) <- rownames(corr.X.Cx) <- Xnoms
+    colnames(corr.Y.Cy) <- colnames(corr.Y.Cx) <- axenames
+    colnames(corr.X.Cy) <- colnames(corr.X.Cx) <- axenames
+
     ## Compute the two redundancy statistics
     RsquareY.X <- simpleRDA2(Y, X)
     RsquareX.Y <- simpleRDA2(X, Y)
@@ -143,12 +168,13 @@
         p.perm <- NA
     }
 
-    out <- list(Pillai=PillaiTrace, EigenValues=EigenValues, CanCorr=K.svd$d,
+    out <- list(Pillai=PillaiTrace, Eigenvalues=Eigenvalues, CanCorr=K.svd$d,
                 Mat.ranks=c(RsquareX.Y$m, RsquareY.X$m), 
                 RDA.Rsquares=c(RsquareY.X$Rsquare, RsquareX.Y$Rsquare),
                 RDA.adj.Rsq=c(Rsquare.adj.Y.X, Rsquare.adj.X.Y),
-                nperm=nperm, p.Pillai=p.Pillai, p.perm=p.perm,
-                AA=AA, BB=BB, Cy=Cy, Cx=Cx, call = match.call())
+                nperm=nperm, p.Pillai=p.Pillai, p.perm=p.perm, Cy=Cy, Cx=Cx, 
+                corr.Y.Cy=corr.Y.Cy, corr.X.Cx=corr.X.Cx, corr.Y.Cx=corr.Y.Cx, 
+                corr.X.Cy=corr.X.Cy, call = match.call())
     class(out) <- "CCorA"
     out
 }

Modified: branches/1.17/R/biplot.CCorA.R
===================================================================
--- branches/1.17/R/biplot.CCorA.R	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/R/biplot.CCorA.R	2010-02-16 06:47:38 UTC (rev 1126)
@@ -1,25 +1,121 @@
 `biplot.CCorA` <-
-    function(x, xlabs, which = 1:2, ...)
+    function(x, plot.type="ov", xlabs, plot.axes = 1:2, int=0.5, col.Y="red", col.X="blue", cex=c(0.7,0.9), ...)
 {
-    plottable <- x$Mat.ranks[which] > 1
-    if (!all(plottable)) {
-        warning("plot ", paste(which(!plottable), collapse=","),
-                " not drawn because it only has one dimension")
-        which <- which[plottable]
-    }
-    if (prod(par("mfrow")) < length(which)) {
-        op <- par(mfrow=c(1,2))
-        on.exit(par(op))
-    }
-    if (missing(xlabs))
-        xlabs <- rownames(x$Cy)
-    else if (!is.null(xlabs) && is.na(xlabs))
-        xlabs <- rep(NA, nrow(x$Cy))
-    else if (is.null(xlabs))
-        xlabs <- 1:nrow(x$Cy)
-    if (any(which == 1))
-        biplot(x$Cy, x$AA, xlabs = xlabs, ...)
-    if (any(which == 2))
-        biplot(x$Cx, x$BB, xlabs = xlabs, ...)
-    invisible()
-}
+    #### Internal function
+	larger.frame <- function(mat, percent=0.10)
+	# Produce an object plot 10% larger than strictly necessary
+	{
+	range.mat <- apply(mat,2,range)
+	z <- apply(range.mat, 2, function(x) x[2]-x[1])
+	range.mat[1,] <- range.mat[1,]-z*percent
+	range.mat[2,] <- range.mat[2,]+z*percent
+	range.mat
+	}
+	####
+	
+	TYPE <- c("objects","variables","ov","biplots")
+	type <- pmatch(plot.type, TYPE)
+	if(is.na(type)) stop("Invalid plot.type")
+
+    epsilon <- sqrt(.Machine$double.eps)
+	if(length(which(x$Eigenvalues > epsilon)) == 1)
+		stop("Plot of axes (", paste(plot.axes, collapse=","),
+			") not drawn because the solution has a single dimension.")
+	if(max(plot.axes) > length(which(x$Eigenvalues > epsilon)))
+		stop("Plot of axes (", paste(plot.axes, collapse=","),
+			") not drawn because the solution has fewer dimensions.")
+
+	if (missing(xlabs))
+		xlabs <- rownames(x$Cy)
+	else if (!is.null(xlabs) && is.na(xlabs))
+		xlabs <- rep(NA, nrow(x$Cy))
+	else if (is.null(xlabs))
+		xlabs <- 1:nrow(x$Cy)
+	#
+	lf.Y <- larger.frame(x$Cy[,plot.axes])
+	lf.X <- larger.frame(x$Cx[,plot.axes])
+	#
+	# Four plot types are available
+	if(type == 1) {  # Object plots
+		# cat('plot.type = objects')
+		par(mfrow=c(1,2), pty = "s")
+		plot(lf.Y, asp=1, xlab=colnames(x$Cy[,plot.axes[1]]), ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		points(x$Cy[,plot.axes], col=col.Y)   # Solid dot: pch=19 
+		text(x$Cy[,plot.axes],labels=xlabs, pos=3, col=col.Y, cex=cex[1])
+		title(main = c("CCorA object plot","First data table (Y)"), line=2)
+		#
+		plot(lf.X, asp=1, xlab=colnames(x$Cy[,plot.axes[1]]), ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		points(x$Cx[,plot.axes], col=col.X)   # Solid dot: pch=19 
+		text(x$Cx[,plot.axes],labels=xlabs, pos=3, col=col.X, cex=cex[1])
+		title(main = c("CCorA object plot","Second data table (X)"), line=2)
+		###
+		###
+		} else if(type == 2) {  # Variable plots
+		# cat('plot.type = variables')
+		par(mfrow=c(1,2), pty = "s")
+		plot(x$corr.Y.Cy[,plot.axes], asp=1, xlim=c(-1,1), ylim=c(-1,1), xlab=colnames(x$Cy[,plot.axes[1]]), 
+			ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		text(x$corr.Y.Cy[,plot.axes],labels=rownames(x$corr.Y.Cy), pos=3, col=col.Y, cex=cex[2])
+		arrows(0,0,x$corr.Y.Cy[,plot.axes[1]],x$corr.Y.Cy[,plot.axes[2]], length=0.05, col=col.Y)
+		abline(h=0, v=0)
+		lines(cos(seq(0, 2*pi, l=100)), sin(seq(0, 2*pi, l=100)))
+		lines(int * cos(seq(0, 2*pi, l=100)), int * sin(seq(0, 2*pi, l=100)))
+		title(main = c("CCorA variable plot","First data table (Y)"), line=2)
+		#
+		plot(x$corr.X.Cx[,plot.axes], asp=1, xlim=c(-1,1), ylim=c(-1,1), xlab=colnames(x$Cy[,plot.axes[1]]), 
+			ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		text(x$corr.X.Cx[,plot.axes],labels=rownames(x$corr.X.Cx), pos=3, col=col.X, cex=cex[2])
+		arrows(0,0,x$corr.X.Cx[,plot.axes[1]],x$corr.X.Cx[,plot.axes[2]], length=0.05, col=col.X)
+		abline(h=0, v=0)
+	    lines(cos(seq(0, 2*pi, l=100)), sin(seq(0, 2*pi, l=100)))
+	    lines(int * cos(seq(0, 2*pi, l=100)), int * sin(seq(0, 2*pi, l=100)))
+		title(main = c("CCorA variable plot","Second data table (X)"), line=2)
+		###
+		###
+		} else if(type == 3) {  # Object and variable plots
+		# cat('plot.type = ov')
+		# par(mfrow=c(2,2), mar=c(4.5,3.5,2,1))
+        layout(matrix(c(1,2,3,4), ncol = 2, nrow = 2, 
+            byrow = TRUE), widths = 1, heights = c(0.5,0.5))
+        par(pty = "s", mar = c(4.5,3.5,2,1))
+		#
+		plot(lf.Y, asp=1, xlab=colnames(x$Cy[,plot.axes[1]]), ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		points(x$Cy[,plot.axes], col=col.Y)   # Solid dot: pch=19 
+		text(x$Cy[,plot.axes],labels=xlabs, pos=3, col=col.Y, cex=cex[1])
+		title(main = c("First data table (Y)"), line=1)
+		#
+		plot(lf.X, asp=1, xlab=colnames(x$Cy[,plot.axes[1]]), ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		points(x$Cx[,plot.axes], col=col.X)   # Solid dot: pch=19 
+		text(x$Cx[,plot.axes],labels=xlabs, pos=3, col=col.X, cex=cex[1])
+		title(main = c("Second data table (X)"), line=1)
+		#
+		plot(x$corr.Y.Cy[,plot.axes], asp=1, xlim=c(-1,1), ylim=c(-1,1), xlab=colnames(x$Cy[,plot.axes[1]]), 
+			ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		text(x$corr.Y.Cy[,plot.axes],labels=rownames(x$corr.Y.Cy), pos=3, col=col.Y, cex=cex[2])
+		arrows(0,0,x$corr.Y.Cy[,plot.axes[1]],x$corr.Y.Cy[,plot.axes[2]], length=0.05, col=col.Y)
+		abline(h=0, v=0)
+		lines(cos(seq(0, 2*pi, l=100)), sin(seq(0, 2*pi, l=100)))
+		lines(int * cos(seq(0, 2*pi, l=100)), int * sin(seq(0, 2*pi, l=100)))
+		#
+		plot(x$corr.X.Cx[,plot.axes], asp=1, xlim=c(-1,1), ylim=c(-1,1), xlab=colnames(x$Cy[,plot.axes[1]]), 
+			ylab=colnames(x$Cy[,plot.axes[2]]), type="n")
+		text(x$corr.X.Cx[,plot.axes],labels=rownames(x$corr.X.Cx), pos=3, col=col.X, cex=cex[2])
+		arrows(0,0,x$corr.X.Cx[,plot.axes[1]],x$corr.X.Cx[,plot.axes[2]], length=0.05, col=col.X)
+		abline(h=0, v=0)
+	    lines(cos(seq(0, 2*pi, l=100)), sin(seq(0, 2*pi, l=100)))
+	    lines(int * cos(seq(0, 2*pi, l=100)), int * sin(seq(0, 2*pi, l=100)))
+		###
+		###
+		} else if(type == 4) {  # Biplots
+		# cat('plot.type = biplot')
+		par(mfrow=c(1,2), pty = "s")
+		biplot(x$Cy[,plot.axes], x$corr.Y.Cy[,plot.axes], col=c("black",col.Y), xlim=lf.Y[,1], ylim=lf.Y[,2], 
+			xlabs = xlabs, arrow.len=0.05, cex=cex, ...)
+		title(main = c("CCorA biplot","First data table (Y)"), line=4)
+		#
+		biplot(x$Cx[,plot.axes], x$corr.X.Cx[,plot.axes], col=c("black",col.X), xlim=lf.X[,1], ylim=lf.X[,2], 
+			xlabs = xlabs, arrow.len=0.05, cex=cex, ...)
+		title(main = c("CCorA biplot","Second data table (X)"), line=4)
+		}
+	invisible()
+}
\ No newline at end of file

Modified: branches/1.17/R/envfit.default.R
===================================================================
--- branches/1.17/R/envfit.default.R	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/R/envfit.default.R	2010-02-16 06:47:38 UTC (rev 1126)
@@ -36,13 +36,15 @@
             if (!is.null(seed)) {
                 assign(".Random.seed", seed, envir = .GlobalEnv)
             }
-            factors <- factorfit(X, Pfac, permutations, strata,  ...)
+            factors <- factorfit(X, Pfac, permutations, strata,
+                                 choices, w = w, ...)
             sol <- list(vector = vectors, factors = factors)
         }
         else vectors <- vectorfit(X, env, permutations, strata, 
                                   choices, w = w, ...)
     }
-    else vectors <- vectorfit(X, env, permutations, strata, ...)
+    else vectors <- vectorfit(X, env, permutations, strata,
+                              choices, w = w, ...)
     sol <- list(vectors = vectors, factors = factors)
     if (!is.null(na.action))
         sol$na.action <- na.action

Modified: branches/1.17/R/print.CCorA.R
===================================================================
--- branches/1.17/R/print.CCorA.R	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/R/print.CCorA.R	2010-02-16 06:47:38 UTC (rev 1126)
@@ -16,7 +16,7 @@
     }
     cat("from F-distribution: ", format.pval(x$p.Pillai), "\n\n")
     out <- rbind("Eigenvalues" = x$EigenValues, "Canonical Correlations" = x$CanCorr)
-    colnames(out) <- colnames(x$AA)
+    colnames(out) <- colnames(x$Cy)
     printCoefmat(out, ...)
     cat("\n")
     out <- rbind("RDA R squares" = x$RDA.Rsquares, "adj. RDA R squares" = x$RDA.adj.Rsq)

Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/inst/ChangeLog	2010-02-16 06:47:38 UTC (rev 1126)
@@ -2,6 +2,13 @@
 
 VEGAN RELEASE VERSIONS at http://cran.r-project.org/
 
+Version 1.17-1 (opened February 16, 2010)
+
+	* merged r1124: envfit() weights fixed for cca().
+
+	* merged r1123: CCorA() output fixes plus upgrades of
+	biplot.CCorA(). 
+	
 Version 1.17-0 (released January 11, 2010)
 
 	* Copied from vegan/pkg/vegan version 1.16-34 (rev1087).


Property changes on: branches/1.17/man
___________________________________________________________________
Name: svn:mergeinfo
   - 

Modified: branches/1.17/man/CCorA.Rd
===================================================================
--- branches/1.17/man/CCorA.Rd	2010-02-15 06:51:19 UTC (rev 1125)
+++ branches/1.17/man/CCorA.Rd	2010-02-16 06:47:38 UTC (rev 1126)
@@ -8,31 +8,48 @@
 
 \description{Canonical correlation analysis, following Brian McArdle's
 unpublished graduate course notes, plus improvements to allow the
-calculations in the case of very sparse and collinear matrices. }
+calculations in the case of very sparse and collinear matrices, and 
+permutation test of Pillai's trace statistic. }
 
 \usage{
 CCorA(Y, X, stand.Y=FALSE, stand.X=FALSE, nperm = 0, ...)
 
-\method{biplot}{CCorA}(x, xlabs, which = 1:2, ...)
+\method{biplot}{CCorA}(x, plot.type="ov", xlabs, plot.axes = 1:2, int=0.5, 
+   col.Y="red", col.X="blue", cex=c(0.7,0.9), ...)
 }
 
 \arguments{
-  \item{Y}{ left matrix. }
-  \item{X}{ right matrix. }
-  \item{stand.Y}{ logical; should \code{Y} be standardized? }
-  \item{stand.X}{ logical; should \code{X} be standardized? }
-  \item{nperm}{ numeric; Number of permutations to evaluate the
-    significance of Pillai's trace}
-  \item{x}{\code{CCoaR} result object}
-  \item{xlabs}{Row labels. The default is to use row names, \code{NULL}
+  \item{Y}{ Left matrix (object class: \code{matrix} or \code{data.frame}). }
+  \item{X}{ Right matrix (object class: \code{matrix} or \code{data.frame}). }
+  \item{stand.Y}{ Logical; should \code{Y} be standardized? }
+  \item{stand.X}{ Logical; should \code{X} be standardized? }
+  \item{nperm}{ Numeric; number of permutations to evaluate the
+    significance of Pillai's trace, e.g. \code{nperm=99} or \code{nperm=999}.}
+  \item{x}{\code{CCoaR} result object.}
+  \item{plot.type}{ A character string indicating which of the following 
+    plots should be produced: \code{"objects"}, \code{"variables"}, \code{"ov"} 
+    (separate graphs for objects and variables), or \code{"biplots"}. Any 
+    unambiguous subset containing the first letters of these names can be used 
+    instead of the full names. }
+  \item{xlabs}{ Row labels. The default is to use row names, \code{NULL}
     uses row numbers instead, and \code{NA} suppresses plotting row names
-    completely}
-  \item{which}{ \code{1} plots \code{Y} results, and
-    \code{2} plots \code{X1} results }
-  \item{\dots}{Other arguments passed to functions. \code{biplot.CCorA}
-    passes graphical arguments to \code{\link{biplot}} and
-    \code{\link{biplot.default}}, \code{CCorA} currently ignores extra 
-    arguments.} 
+    completely.}
+  \item{plot.axes}{ A vector with 2 values containing the order numbers of 
+    the canonical axes to be plotted. Default: first two axes. }
+  \item{int}{ Radius of the inner circles plotted as visual references in 
+    the plots of the variables. Default: \code{int=0.5}. With \code{int=0}, 
+    no inner circle is plotted. }
+  \item{col.Y}{ Color used for objects and variables in the first data 
+    table (Y) plots. In biplots, the objects are in black. }
+  \item{col.X}{ Color used for objects and variables in the second data 
+    table (X) plots. }
+  \item{cex}{ A vector with 2 values containing the size reduction factors 
+    for the object and variable names, respectively, in the plots. 
+    Default values: \code{cex=c(0.7,0.9)}. }
+  \item{\dots}{ Other arguments passed to these functions. The function 
+    \code{biplot.CCorA} passes graphical arguments to \code{\link{biplot}} 
+    and \code{\link{biplot.default}}. \code{CCorA} currently ignores extra 
+    arguments. } 
 }
 
 \details{
@@ -40,58 +57,61 @@
   combinations of the variables of \code{Y} that are maximally
   correlated to linear combinations of the variables of \code{X}. The
   analysis estimates the relationships and displays them in graphs.
+  Pillai's trace statistic is computed and tested parametrically (F-test);
+  a permutation test is also available.
 
-  Algorithmic notes:
-  \enumerate{
-    \item
-    All data matrices are replaced by their PCA object scores, computed
-    by SVD.
-    \item
-    The blunt approach would be to read the three matrices, compute the
-    covariance matrices, then the matrix
-    \code{S12 \%*\% inv(S22) \%*\% t(S12) \%*\% inv(S11)}.
-    Its trace is Pillai's trace statistic. 
-    \item
-    This approach may fail, however, when there is heavy multicollinearity
-    in very sparse data matrices, as it is the case in 4th-corner inflated
-    data matrices for example. The safe approach is to replace all data
-    matrices by their PCA object scores.
-    \item
-    Inversion by \code{\link{solve}} is avoided. Computation of inverses
-    is done by \acronym{SVD}  (\code{\link{svd}}) in most cases.
-    \item
-    Regression by \acronym{OLS} is also avoided. Regression residuals are
-    computed by \acronym{QR} decomposition (\code{\link{qr}}).
-  }
+  Algorithmic note -- 
+  The blunt approach would be to read the two matrices, compute the
+  covariance matrices, then the matrix
+  \code{S12 \%*\% inv(S22) \%*\% t(S12) \%*\% inv(S11)}.
+  Its trace is Pillai's trace statistic. 
+  This approach may fail, however, when there is heavy multicollinearity
+  in very sparse data matrices. The safe approach is to replace all data
+  matrices by their PCA object scores.
 
-The \code{biplot} function can produce two biplots, each for the left
-matrix and right matrix solutions. The function passes all arguments to
-\code{\link{biplot.default}}, and you should consult its help page for
-configuring biplots.
+The function can produce different types of plots depending on the option 
+chosen: 
+\code{"objects"} produces two plots of the objects, one in the space 
+of Y, the second in the space of X; 
+\code{"variables"} produces two plots of the variables, one of the variables 
+of Y in the space of Y, the second of the variables of X in the space of X; 
+\code{"ov"} produces four plots, two of the objects and two of the variables; 
+\code{"biplots"} produces two biplots, one for the first matrix (Y) and 
+one for second matrix (X) solutions. For biplots, the function passes all arguments 
+to \code{\link{biplot.default}}; consult its help page for configuring biplots.
 }
 
 \value{
-  Function \code{CCorA} returns a list containing the following components:
-  \item{ Pillai }{ Pillai's trace statistic = sum of canonical
+  Function \code{CCorA} returns a list containing the following elements:
+  \item{ Pillai }{ Pillai's trace statistic = sum of the canonical
     eigenvalues. } 
-  \item{ EigenValues }{ Canonical eigenvalues. They are the squares of the
+  \item{ Eigenvalues }{ Canonical eigenvalues. They are the squares of the
     canonical correlations. }
   \item{ CanCorr }{ Canonical correlations. }
-  \item{ Mat.ranks }{ Ranks of matrices \code{Y} and X1 (possibly after
-    controlling for X2). }
+  \item{ Mat.ranks }{ Ranks of matrices \code{Y} and \code{X}. }
   \item{ RDA.Rsquares }{ Bimultivariate redundancy coefficients
-    (R-squares) of RDAs of Y|X1 and X1|Y. } 
-  \item{ RDA.adj.Rsq }{ \code{RDA.Rsquares} adjusted for n and number of
-    explanatory variables. }
-  \item{ AA }{ Scores of Y variables in Y biplot. }
-  \item{ BB }{ Scores of X1 variables in X1 biplot. }
+    (R-squares) of RDAs of Y|X and X|Y. } 
+  \item{ RDA.adj.Rsq }{ \code{RDA.Rsquares} adjusted for \code{n} and the number 
+    of explanatory variables. }
+  \item{ nperm }{ Number of permutations. }
+  \item{ p.Pillai }{ Parametric probability value associated with Pillai's trace. }
+  \item{ p.perm }{ Permutational probability associated with Pillai's trace. }
   \item{ Cy }{ Object scores in Y biplot. }
-  \item{ Cx }{ Object scores in X1 biplot. }
+  \item{ Cx }{ Object scores in X biplot. }
+  \item{ corr.Y.Cy }{ Scores of Y variables in Y biplot, computed as cor(Y,Cy). }
+  \item{ corr.X.Cx }{ Scores of X variables in X biplot, computed as cor(X,Cx). }
+  \item{ corr.Y.Cx }{ cor(Y,Cy) available for plotting variables Y in space of X manually. }
+  \item{ corr.X.Cy }{ cor(X,Cx) available for plotting variables X in space of Y manually. }
+  \item{ call }{ Call to the CCorA function. }
 }
 
 \references{ 
   Hotelling, H. 1936. Relations between two sets of
   variates. \emph{Biometrika} \strong{28}: 321-377.
+  
+  Legendre, P. 2005. Species associations: the Kendall coefficient of 
+  concordance revisited. \emph{Journal of Agricultural, Biological, and 
+  Environmental Statistics} \strong{10}: 226-245.
 }
 
 \author{ Pierre Legendre, Departement de Sciences Biologiques,
@@ -99,18 +119,32 @@
   Jari Oksanen. }
 
 \examples{
-# Example using random numbers
+# Example using two mite groups. The mite data are available in vegan
+data(mite)
+# Two mite species associations (Legendre 2005, Fig. 4)
+group.1 <- c(1,2,4:8,10:15,17,19:22,24,26:30)
+group.2 <- c(3,9,16,18,23,25,31:35)
+# Separate Hellinger transformations of the two groups of species 
+mite.hel.1 <- decostand(mite[,group.1], "hel")
+mite.hel.2 <- decostand(mite[,group.2], "hel")
+rownames(mite.hel.1) = paste("S",1:nrow(mite),sep="")
+rownames(mite.hel.2) = paste("S",1:nrow(mite),sep="")
+out <- CCorA(mite.hel.1, mite.hel.2)
+out
+summary(out)
+biplot(out, "ob")                 # Two plots of objects
+biplot(out, "v", cex=c(0.7,0.6))  # Two plots of variables
+biplot(out, "ov", cex=c(0.7,0.6)) # Four plots (2 for objects, 2 for variables)
+biplot(out, "b", cex=c(0.7,0.6))  # Two biplots
+biplot(out, xlabs = NA, plot.axes = c(3,5))    # Plot axes 3, 5. No object names
+biplot(out, plot.type="biplots", xlabs = NULL) # Replace object names by numbers
+
+# Example using random numbers. No significant relationship is expected
 mat1 <- matrix(rnorm(60),20,3)
 mat2 <- matrix(rnorm(100),20,5)
-CCorA(mat1, mat2)
-
-# Example using intercountry life-cycle savings data, 50 countries
-data(LifeCycleSavings)
-pop <- LifeCycleSavings[, 2:3]
-oec <- LifeCycleSavings[, -(2:3)]
-out <- CCorA(pop, oec)
-out
-biplot(out, xlabs = NA)
+out2 = CCorA(mat1, mat2, nperm=99)
+out2
+biplot(out2, "b")
 }
 
 \keyword{ multivariate }



More information about the Vegan-commits mailing list