R/biplot.pcoa.R

   1 'biplot.pcoa' <-
   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
   9 #
  10 # Author: Pierre Legendre, January 2009
  11         {
  12         k <- ncol(x$vectors)
  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.")
  16
  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
  21
  22         if(is.null(Y)) {
  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))
  28
  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)
  33
  34                 } else {
  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)
  38                 n <- nrow(Y)
  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])
  43
  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)
  47         }
  48     invisible()
  49 }