[R] reflecting a PCA biplot

[Gavin Simpson]

On Tue, 2011-08-09 at 22:57 +1000, Andrew Halford wrote:
> Hi Listers,
> 
> I am trying to reflect a PCA biplot in the x-axis (i.e. PC1) but am not
> having much success. In theory I believe all I need to do is multiply the
> site and species scores for the PC1 by -1, which would effectively flip the
> biplot.
> 
> I am creating a blank plot using the plot command and accessing the results
> from a call to rda. I then use the calls to scores to obtain separate site
> and species coordinates and I have worked out how to multiply the
> appropriate PC1 scores by -1 to create the site and species scores I want.
> However I am not sure how to change the call to plot which accesses the
> results of the call to rda to draw the blank plot. The coordinates it is
> accessing are for the unreflected ordination and this does not match the new
> site and species scores that I have.
> 
> 
> > fish.pca <- rda(fish.hel)
> > fish.site <- scores(fish.pca,display="sites",scaling=3)
> > fish.spp <- scores(fish.pca,display="species",scaling=3)
> 
> > fish.site[,"PC1"] <- -1*(fish.site[,"PC1"])
> > fish.spp[,"PC1"] <- -1*(fish.spp[,"PC1"])
> 
> > graph <- plot(fish.pca,display=c("sites","species"),type="n",scaling=3) #
> how do I get the plot to draw up the blank display based on the reversed
> site and species scores?

Do you mean something like...?

require(vegan)
data(dune)
mod <- rda(dune)
si.scrs <- scores(mod, display = "sites", scaling = 3, choices = 1:2)
sp.scrs <- scores(mod, display = "species", scaling = 3, choices = 1:2)

si.scrs[, "PC1"] <- -1 * si.scrs[, "PC1"]
sp.scrs[, "PC1"] <- -1 * sp.scrs[, "PC1"]

xlim <- range(0, si.scrs[, 1], sp.scrs[, 1])
ylim <- range(0, si.scrs[, 2], sp.scrs[, 2])

plot(si.scrs[,1], si.scrs[,2], ylim = ylim, xlim = xlim,
     ylab = "PC2", xlab = "PC1", cex = 0.7, asp = 1)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(sp.scrs[,1], sp.scrs[,2], col = "red", pch = 3, cex = 0.7)
box()

For non-standard plotting of ordination objects, our advice has always
been to build the plot up from lower-level plotting functions rather
than the specific methods supplied with vegan.

A comparison: (not quite the same, I know, but close enough)

layout(matrix(1:2, ncol = 2))
plot(mod, display=c("sites","species"), type = "p", scaling=3,
     main = "Original")
plot(si.scrs[,1], si.scrs[,2], ylim = ylim, xlim = xlim,
     ylab = "PC2", xlab = "PC1", cex = 0.7, asp = 1,
     main = "Flipped PC1")
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(sp.scrs[,1], sp.scrs[,2], col = "red", pch = 3, cex = 0.7)
box()
layout(1)

If you want type = "t", the default for plot.cca, then use the same
coordinates but extract the relevant labels from the two scores objects:

plot(si.scrs[,1], si.scrs[,2], ylim = ylim, xlim = xlim,
     ylab = "PC2", xlab = "PC1", cex = 0.7, asp = 1,
     type = "n") ## no plotting this time first
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
text(si.scrs[,1], si.scrs[,2], labels = rownames(si.scrs),
     cex = 0.7) ## add site and species scores as labels
text(sp.scrs[,1], sp.scrs[,2], labels = rownames(sp.scrs),
     col = "red", cex = 0.7)
box()

HTH

G

> Any help appreciated.
> 
> cheers
> 
> Andy

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%