2 function(x, Y=NULL, plot.axes = c(1,2), dir.axis1=1, dir.axis2=1, rn=NULL, ...)
3 # x = output object from function pcoa.R
4 # Y = optional sites-by-variables data table
5 # plot.axes = the two axes to be plotted
6 # rn = an optional vector, length n, of object labels
7 # dir.axis.1 = -1 to revert axis 1 for the projection of points and variables
8 # dir.axis.2 = -1 to revert axis 2 for the projection of points and variables
10 # Author: Pierre Legendre, January 2009
13 if(k < 2) stop("There is a single eigenvalue. No plot can be produced.")
14 if(k < plot.axes[1]) stop("Axis",plot.axes[1],"does not exist.")
15 if(k < plot.axes[2]) stop("Axis",plot.axes[2],"does not exist.")
17 if(!is.null(rn)) rownames(x$vectors) <- rn
18 labels = colnames(x$vectors[,plot.axes])
19 diag.dir <- diag(c(dir.axis1,dir.axis2))
20 x$vectors[,plot.axes] <- x$vectors[,plot.axes] %*% diag.dir
23 limits <- apply(x$vectors[,plot.axes], 2, range)
24 ran.x <- limits[2,1] - limits[1,1]
25 ran.y <- limits[2,2] - limits[1,2]
26 xlim <- c((limits[1,1]-ran.x/10), (limits[2,1]+ran.x/5))
27 ylim <- c((limits[1,2]-ran.y/10), (limits[2,2]+ran.y/10))
29 par(mai = c(1.0, 1.0, 1.0, 0.5))
30 plot(x$vectors[,plot.axes], xlab=labels[1], ylab=labels[2], xlim=xlim, ylim=ylim, asp=1)
31 text(x$vectors[,plot.axes], labels=rownames(x$vectors), pos=4, cex=1, offset=0.5)
32 title(main = "PCoA ordination", line=2.5)
35 # Find positions of variables in biplot:
36 # construct U from covariance matrix between Y and standardized point vectors
37 # (equivalent to PCA scaling 1, since PCoA preserves distances among objects)
39 points.stand <- scale(x$vectors[,plot.axes])
40 S <- cov(Y, points.stand)
41 U <- S %*% diag((x$values$Eigenvalues[plot.axes]/(n-1))^(-0.5))
42 colnames(U) <- colnames(x$vectors[,plot.axes])
44 par(mai = c(1, 0.5, 1.4, 0))
45 biplot(x$vectors[,plot.axes], U, xlab=labels[1], ylab=labels[2])
46 title(main = c("PCoA biplot","Response variables projected","as in PCA with scaling 1"), line=4)