[Vegan-commits] r1123 - in pkg/vegan: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 14 08:41:49 CET 2010
Author: jarioksa
Date: 2010-02-14 08:41:48 +0100 (Sun, 14 Feb 2010)
New Revision: 1123
Modified:
pkg/vegan/R/CCorA.R
pkg/vegan/R/biplot.CCorA.R
pkg/vegan/R/print.CCorA.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/CCorA.Rd
Log:
CCorA: fixes bugs in presentation of results and adds new biplot types (P Legendre)
Modified: pkg/vegan/R/CCorA.R
===================================================================
--- pkg/vegan/R/CCorA.R 2010-02-10 17:36:16 UTC (rev 1122)
+++ pkg/vegan/R/CCorA.R 2010-02-14 07:41:48 UTC (rev 1123)
@@ -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: pkg/vegan/R/biplot.CCorA.R
===================================================================
--- pkg/vegan/R/biplot.CCorA.R 2010-02-10 17:36:16 UTC (rev 1122)
+++ pkg/vegan/R/biplot.CCorA.R 2010-02-14 07:41:48 UTC (rev 1123)
@@ -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: pkg/vegan/R/print.CCorA.R
===================================================================
--- pkg/vegan/R/print.CCorA.R 2010-02-10 17:36:16 UTC (rev 1122)
+++ pkg/vegan/R/print.CCorA.R 2010-02-14 07:41:48 UTC (rev 1123)
@@ -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: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2010-02-10 17:36:16 UTC (rev 1122)
+++ pkg/vegan/inst/ChangeLog 2010-02-14 07:41:48 UTC (rev 1123)
@@ -4,6 +4,10 @@
Version 1.18-1 (opened February 9, 2010)
+ * CCorA: Fixed bug in presentation of variables in plots. Adds new
+ biplot types. General improvement in checking exceptional cases
+ improve stability.
+
* predict.cca, predict.rda: match 'newdata' by row names or column
names in type = "wa" and type = "sp". This is similar as
predict.prcomp/princomp. Gained choice type = "working" for
Modified: pkg/vegan/man/CCorA.Rd
===================================================================
--- pkg/vegan/man/CCorA.Rd 2010-02-10 17:36:16 UTC (rev 1122)
+++ pkg/vegan/man/CCorA.Rd 2010-02-14 07:41:48 UTC (rev 1123)
@@ -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