Reputation: 679
I'm trying to produce pca biplots for data subsets. Within the same principal components environment I'd like to plot only subsets based on Moisture levels.
# Packages
library(vegan)
# Sample data
data(dune, dune.env)
dd <- cbind(dune.env, dune)
# Runnig PCA
pca1 <- rda(dd[, -c(1:5)], scale=T)
# Plot
plot(pca1, scaling=3)
# Now, instead the plot above, I'd like to get 5 different plots (one per dd$Moisture level) where I see the principal components scores but only for the subsets based on:
levels(dd$Moisture)
Many thanks in advance for any inputs!!
Upvotes: 2
Views: 2190
Reputation: 174813
There are somewhat easier ways to do these faceted plots with ordixyplot()
in vegan or the ggvegan package (currently alpha so simple things sometimes need to be done by hand).
# Packages
library("vegan")
library("ggvegan")
# Sample data
data(dune, dune.env)
# PCA
pca1 <- rda(dune, scale = TRUE)
As I said, you need to do a bit of munging by hand and the moment, but you get the legend for free with ggplot.
## use fortify method to extract scores in ggplot-friendly format
scrs <- fortify(pca1, scaling = 3)
## take only site scores for this
scrs <- with(scrs, scrs[Score == "sites", ])
## add on Moisture variable
scrs <- cbind(scrs, Moisture = dune.env$Moisture)
## for all points on one plot, Moisture coded by colour
plt <- ggplot(scrs, aes(x = Dim1, y = Dim2, colour = Moisture)) +
geom_point() + coord_fixed()
plt
## to facet that plot into multiple panels
plt + facet_wrap( ~ Moisture)
ordixyplot()
versionMore of the work is canned within ordixyplot()
than with the ggvegan version, but you have to work a little harder to add the key (legend) and I can never remember how to do this with Lattice.
## Single plot
ordixyplot(pca1, data = dune.env, formula = PC1 ~ PC2, group = Moisture)
## Facet plot
ordixyplot(pca1, data = dune.env, formula = PC1 ~ PC2 | Moisture)
For base graphics there is an easier version to colour the points on a single plot. One version is
scrs <- scores(pca1, display = "sites")
cols <- c("red","green","blue","orange","purple")
plot(scrs[,1], scrs[,2], pch = 19,
col = cols[as.numeric(as.character(dune.env$Moisture))])
legend("topright", legend = 1:5, title = "Moisture", pch = 19,
col = cols, bty = "n")
You can find out more about doing the plots this way using base graphics in a blog post that I wrote a few years ago.
Upvotes: 1
Reputation: 679
# Packages
library("vegan")
# Preparing data
data(dune, dune.env)
dd <- cbind(dune.env, dune)
# Order data frame based on dd$Moisture
dd <- dd[order(dd$Moisture),]
str(dd)
# Runnig PCA
pca1 <- rda(dd[, -c(1:5)], scale=T)
# Plot
biplot(pca1, scaling=3, type = c("text", "points"))
# I need to get 5 different plots (one per dd$Moisture level)
# grab the scores on axes required
site.scr <- scores(pca1, display = "sites", scaling = 3) # x
spp.scr <- scores(pca1, display = "species", scaling = 3) # y
## generate x,y lims
xlims <- range(site.scr[,1], spp.scr[,1])
ylims <- range(site.scr[,2], spp.scr[,2])
## plot what you want, but leave out sites
biplot(mod, display = "species",
xlim=xlims, ylim=ylims,
scaling = 3)
## now add in scores as you want, to simulate we plot:
# Moisture - All together but can be in independetn plots
points(site.scr[1:7,], col = "blue", pch = 2)
points(site.scr[8:11,], col = "green", pch = 2)
points(site.scr[0:0,], col = "orange", pch = 2) # Level 3 not present
points(site.scr[12:13,], col = "grey", pch = 2)
points(site.scr[14:20,], col = "black", pch = 2)
Upvotes: 0