# 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