# Examples use residuals from a regression of salamander 
# morphological traits against body size (snout to vent length, SVL).
# Observations are species means and a phylogenetic covariance matrix
# describes the relatedness among observations.
data("PlethMorph")
Y <- as.data.frame(PlethMorph[c("TailLength", "HeadLength", 
"Snout.eye", "BodyWidth", 
"Forelimb", "Hindlimb")])
Y <- as.matrix(Y)
R <- lm.rrpp(Y ~ SVL, data = PlethMorph, 
iter = 0, print.progress = FALSE)$LM$residuals
# PCA (on correlation matrix)
PCA.ols <- ordinate(R, scale. = TRUE)
PCA.ols$rot
prcomp(R, scale. = TRUE)$rotation # should be the same
# phyPCA (sensu Revell, 2009)
# with projection of untransformed residuals (Collyer & Adams 2020)
PCA.gls <- ordinate(R, scale. = TRUE, 
transform. = FALSE, 
Cov = PlethMorph$PhyCov)
# phyPCA with transformed residuals (orthogonal projection, 
# Collyer & Adams 2020)
PCA.t.gls <- ordinate(R, scale. = TRUE, 
transform. = TRUE, 
Cov = PlethMorph$PhyCov)
 
 # Align to phylogenetic signal (in each case)
 
 PaCA.ols <- ordinate(R, A = PlethMorph$PhyCov, scale. = TRUE)
 
 PaCA.gls <- ordinate(R, A = PlethMorph$PhyCov, 
 scale. = TRUE,
 transform. = FALSE, 
 Cov = PlethMorph$PhyCov)
 
 PaCA.t.gls <- ordinate(R, A = PlethMorph$PhyCov, 
 scale. = TRUE,
 transform. = TRUE, 
 Cov = PlethMorph$PhyCov)
 
 # Summaries
 
 summary(PCA.ols)
 summary(PCA.gls)
 summary(PCA.t.gls)
 summary(PaCA.ols)
 summary(PaCA.gls)
 summary(PaCA.t.gls)
 
 # Plots
 
 par(mfrow = c(2,3))
 plot(PCA.ols, main = "PCA OLS")
 plot(PCA.gls, main = "PCA GLS")
 plot(PCA.t.gls, main = "PCA t-GLS")
 plot(PaCA.ols, main = "PaCA OLS")
 plot(PaCA.gls, main = "PaCA GLS")
 plot(PaCA.t.gls, main = "PaCA t-GLS")
 par(mfrow = c(1,1))
 
 # Changing some plot aesthetics (the arguments in plot.ordinate and 
 # plot.default are important for changing plot parameters)
 
 P1 <- plot(PaCA.gls, main = "PaCA GLS", include.axes = TRUE)
 
 P2 <- plot(PaCA.gls, main = "PaCA GLS", include.axes = TRUE, 
 col = 4, pch = 21, bg = PlethMorph$group)
 add.tree(P2, PlethMorph$tree, edge.col = 4)
 P3 <- plot(PaCA.gls, main = "PaCA GLS", include.axes = TRUE, 
 col = 4, pch = 21, bg = PlethMorph$group,
 flip = 1)
 add.tree(P3, PlethMorph$tree, edge.col = 4)
 
 
Run the code above in your browser using DataLab