# NOT RUN {
data(plethspecies)
Y.gpa <- gpagen(plethspecies$land) #GPA-alignment
### Traditional PCA
PCA <- gm.prcomp(Y.gpa$coords)
summary(PCA)
plot(PCA, main = "PCA")
plot(PCA, main = "PCA", flip = 1) # flip the first axis
plot(PCA, main = "PCA", axis1 = 3, axis2 = 4) # change PCs viewed
### Phylomorphospace - PCA with phylogeny (result is same as above,
### but with estimated ancestral states projected into plot)
PCA.w.phylo <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy)
summary(PCA.w.phylo)
plot(PCA.w.phylo, phylo = TRUE, main = "PCA.w.phylo")
### Phylogenetic PCA - PCA based on GLS-centering and projection
# This is the same as the method described by Revell (2009)
phylo.PCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, GLS = TRUE)
summary(phylo.PCA)
plot(phylo.PCA, phylo = TRUE, main = "phylo PCA")
### Phylogenetic PCA - PCA based on GLS-centering and transformed
# projection
# This produces a PCA independent of phylogeny
phylo.tPCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy,
GLS = TRUE, transform = TRUE)
summary(phylo.tPCA)
plot(phylo.tPCA, phylo = TRUE, main = "phylo PCA")
### PaCA - Alignment of data to physlogenetic signal rather than axis of
### greatest variation, like in PCA
# OLS method (rotation of PCA)
PaCA.ols <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy,
align.to.phy = TRUE)
summary(PaCA.ols)
plot(PaCA.ols, phylo = TRUE, main = "PaCA using OLS")
# GLS method (rotation of Phylogenetic PCA)
PaCA.gls <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy,
align.to.phy = TRUE, GLS = TRUE)
summary(PaCA.gls)
plot(PaCA.gls, phylo = TRUE, main = "PaCA using GLS")
# GLS method (rotation of Phylogenetic PCA with transformed data)
PaCA.gls <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy,
align.to.phy = TRUE, GLS = TRUE, transform = TRUE)
summary(PaCA.gls)
plot(PaCA.gls, phylo = TRUE,
main = "PaCA using GLS and transformed projection")
### Advanced Plotting
gps <- as.factor(c(rep("gp1", 5), rep("gp2", 4))) # Two random groups
par(mar=c(2, 2, 2, 2))
plot(PaCA.ols, pch=22, cex = 1.5, bg = gps, phylo = TRUE)
# Modify options as desired
# Add things as desired using standard R plotting
text(par()$usr[1], 0.1*par()$usr[3], labels = "PC1 - 45.64%",
pos = 4, font = 2)
text(0, 0.95*par()$usr[4], labels = "PC2 - 18.80%", pos = 4, font = 2)
legend("topleft", pch=22, pt.bg = unique(gps), legend = levels(gps))
### 3D plot with a phylogeny and time on the z-axis
plot(PCA.w.phylo, time.plot = TRUE)
plot(PCA.w.phylo, time.plot = TRUE, bg = "red",
phylo.par = list(tip.labels = TRUE,
tip.txt.cex = 2, edge.color = "blue", edge.width = 2))
# }
Run the code above in your browser using DataLab