Day-02: Principal Components Analysis (PCA) and Visualizing Shape Differences
Content Description
Here we demonstrate some of the graphical capabilities of
geomorph for visualizing patterns of shape variation and
interpreting particular shape differences.
First, let´s load some example data and plot the specimens (before and after GPA alignment):
library(geomorph)
data(plethodon)
Y.gpa <- gpagen(plethodon$land, print.progress = F) # GPA-alignment
par(mfrow=c(1,2))
plotAllSpecimens(plethodon$land, links=plethodon$links) # Raw data
mtext("Raw Data")
plotAllSpecimens(Y.gpa$coords, links=plethodon$links) # GPA-aligned data
mtext("GPA-Aligned Specimens")
par(mfrow=c(1,1)) 1: Principal Components Analysis
The function gm.prcomp provides visualization of
patterns of data in shape space via PCA (and as we’ll see later, via
other data alignment procedures). Here is a simple example of use of the
function:
PCA <- gm.prcomp(Y.gpa$coords)
summary(PCA)##
## Ordination type: Principal Component Analysis
## Centering by OLS mean
## Orthogonal projection of OLS residuals
## Number of observations: 40
## Number of vectors 20
##
## Importance of Components:
## Comp1 Comp2 Comp3 Comp4
## Eigenvalues 0.001855441 0.001566597 0.0004141056 0.0002278444
## Proportion of Variance 0.367433295 0.310233333 0.0820053761 0.0451200584
## Cumulative Proportion 0.367433295 0.677666628 0.7596720044 0.8047920628
## Comp5 Comp6 Comp7 Comp8
## Eigenvalues 0.0001725997 0.0001672575 0.0001549008 0.0001251132
## Proportion of Variance 0.0341799312 0.0331220169 0.0306750315 0.0247761773
## Cumulative Proportion 0.8389719940 0.8720940109 0.9027690424 0.9275452197
## Comp9 Comp10 Comp11 Comp12
## Eigenvalues 8.468503e-05 6.869926e-05 5.430853e-05 4.423148e-05
## Proportion of Variance 1.677019e-02 1.360452e-02 1.075473e-02 8.759165e-03
## Cumulative Proportion 9.443154e-01 9.579199e-01 9.686747e-01 9.774338e-01
## Comp13 Comp14 Comp15 Comp16
## Eigenvalues 0.0000261742 0.0000192644 1.741769e-05 1.715141e-05
## Proportion of Variance 0.0051832795 0.0038149309 3.449226e-03 3.396495e-03
## Cumulative Proportion 0.9826170994 0.9864320303 9.898813e-01 9.932778e-01
## Comp17 Comp18 Comp19 Comp20
## Eigenvalues 1.555244e-05 9.253679e-06 5.336746e-06 3.802714e-06
## Proportion of Variance 3.079853e-03 1.832507e-03 1.056837e-03 7.530520e-04
## Cumulative Proportion 9.963576e-01 9.981901e-01 9.992469e-01 1.000000e+00plot(PCA)
Here we have summarized variation along each PC axis, and provided a plot. Here is a more advanced version:
gps <- as.factor(paste(plethodon$species, plethodon$site)) #define some groups for plotting
plot(PCA, pch=22, cex = 1.5, bg = gps)
legend("topleft", pch=22, pt.bg = unique(gps), legend = levels(gps))
2: Visualize Shape Differences Between Specimens
Taking the previous example, we may now with to generate graphical
depictions of shape differences for certain specimens.
Generating graphical depictions of shapes is accomlished using
plotRefToTarget. Here are a few examples of its use:
Deformations Specimen to Specimen
M <- mshape(Y.gpa$coords)
par(mfrow=c(1,2))
plotRefToTarget(M,Y.gpa$coords[,,39], links=plethodon$links)
mtext("TPS")
plotRefToTarget(M,Y.gpa$coords[,,39], mag=2.5, links=plethodon$links)
mtext("TPS: 2.5X magnification")
par(mfrow=c(1,1))
par(mfrow=c(1,2))
plotRefToTarget(M,Y.gpa$coords[,,39], links=plethodon$links, method="vector", mag=3)
mtext("Vector Displacements")
plotRefToTarget(M,Y.gpa$coords[,,39], links=plethodon$links,gridPars=gridPar(pt.bg="red", link.col="green", pt.size = 1), method="vector", mag=3)
mtext("Vector Displacements: Other Options")
par(mfrow=c(1,1))
par(mfrow=c(1,2))
plotRefToTarget(M,Y.gpa$coords[,,39], mag=2, outline=plethodon$outline)
mtext("Outline Deformation")
plotRefToTarget(M,Y.gpa$coords[,,39], method="points", outline=plethodon$outline)
mtext("Outline Deformations Ref (gray) & and Tar (black)")
par(mfrow=c(1,1))3: Advanced Shape Plotting: Shape Prediction
NOw let’s visualize how shape changes along PC1 of the previous
ordination. For this we must PREDICT shapes
along PC1, and then plot them. This is accomplished using a combination
of the functions shape.predictor (to identify shapes along
PC1), and plotRefToTarget (to visualize them):
PCA Shape Predictions
PC <- PCA$x[,1]
preds <- shape.predictor(Y.gpa$coords, x= PC, Intercept = FALSE,
pred1 = min(PC), pred2 = max(PC)) # PC 1 extremes, more technically
par(mfrow=c(1,2))
plotRefToTarget(M, preds$pred1, links = plethodon$links)
mtext("PC1 - Min.")
plotRefToTarget(M, preds$pred2, links = plethodon$links)
mtext("PC1 - Max.")
par(mfrow=c(1,1))Shape Predictions Along a Regression
Many times, we fit a linear model in the form of regression (e.g.,
allometry). Here is how we visualize shape changes along the extremes of
that regression using shape.predictor:
gdf <- geomorph.data.frame(Y.gpa)
plethAllometry <- procD.lm(coords ~ log(Csize), data=gdf, print.progress = FALSE)
allom.plot <- plot(plethAllometry,
type = "regression",
predictor = log(gdf$Csize),
reg.type ="PredLine") # make sure to have a predictor 
preds <- shape.predictor(plethAllometry$GM$fitted, x= allom.plot$RegScore, Intercept = TRUE,
predmin = min(allom.plot$RegScore),
predmax = max(allom.plot$RegScore))
par(mfrow=c(1,2))
plotRefToTarget(M, preds$predmin, mag=3, links = plethodon$links)
mtext("Regression Min: 3X")
plotRefToTarget(M, preds$predmax, mag=3, links = plethodon$links)
mtext("Regression Max: 3X")
par(mfrow=c(1,1))Shape Predictions of Group Differences
Another very common statistical design is one of evaluating group differences in shape (i.e., ANOVA). Here we wish to visualize shape differences among groups. Below is an example:
gdf <- geomorph.data.frame(Y.gpa, species = plethodon$species, site = plethodon$site)
pleth.anova <- procD.lm(coords ~ species*site, data=gdf, print.progress = FALSE)
X <- pleth.anova$X
X # includes intercept; remove for better functioning ## (Intercept) speciesTeyah siteSymp speciesTeyah:siteSymp
## 1 1 0 1 0
## 2 1 0 1 0
## 3 1 0 1 0
## 4 1 0 1 0
## 5 1 0 1 0
## 6 1 0 1 0
## 7 1 0 1 0
## 8 1 0 1 0
## 9 1 0 1 0
## 10 1 0 1 0
## 11 1 1 1 1
## 12 1 1 1 1
## 13 1 1 1 1
## 14 1 1 1 1
## 15 1 1 1 1
## 16 1 1 1 1
## 17 1 1 1 1
## 18 1 1 1 1
## 19 1 1 1 1
## 20 1 1 1 1
## 21 1 0 0 0
## 22 1 0 0 0
## 23 1 0 0 0
## 24 1 0 0 0
## 25 1 0 0 0
## 26 1 0 0 0
## 27 1 0 0 0
## 28 1 0 0 0
## 29 1 0 0 0
## 30 1 0 0 0
## 31 1 1 0 0
## 32 1 1 0 0
## 33 1 1 0 0
## 34 1 1 0 0
## 35 1 1 0 0
## 36 1 1 0 0
## 37 1 1 0 0
## 38 1 1 0 0
## 39 1 1 0 0
## 40 1 1 0 0X <- X[,-1]
symJord <- c(0,1,0) # design for P. Jordani in sympatry
alloJord <- c(0,0,0) # design for P. Jordani in allopatry
preds <- shape.predictor(arrayspecs(pleth.anova$fitted, 12, 2), x = X, Intercept = TRUE,
symJord=symJord, alloJord=alloJord)
par(mfrow=c(1,2))
plotRefToTarget(M, preds$symJord, links = plethodon$links, mag=2)
mtext("Sympatric P. Jordani: 2X")
plotRefToTarget(M, preds$alloJord, links = plethodon$links, mag=2)
mtext("Allopatric P. Jordani: 2X")
par(mfrow=c(1,1))4: Selecting Shapes to Visualize in Real Time
One recent enhancement of geomorph is the
picknplot.shape function. With this function, one can
select locations in a statistical dataspace, and geomorph
will mathematially back-translate the statistical summary values in the
plot to obtain landmark coordinates, and will then generate a thin-plate
spline deformation grid of a specimen that would exist at that location
in the dataspace.
This real-time visualization is especially useful for viewing regions of morphospace that are unoccupied by the sample of specimens, and for exploring other non-sampled shapes.
Note: because this is an interactive function, the code is repeated below but not run.
# 2d
# data(plethodon)
# Y.gpa <- gpagen(plethodon$land)
# pleth.pca <- gm.prcomp(Y.gpa$coords)
# pleth.pca.plot <- plot(pleth.pca)
# picknplot.shape(pleth.pca.plot)
# May change arguments for plotRefToTarget
# picknplot.shape(plot(pleth.pca), method = "points", mag = 3, links=plethodon$links)5: 3D warping
Geomorph also has the ability to generate shape
deformations of 3D objects. Again this is interactive code, so it is
repeated, but not run, below:
#scallops <- readland.tps("Data/scallops for viz.tps", specID = "ID")
#ref <- mshape(scallops)
#refmesh <- warpRefMesh(read.ply("Data/glyp02L.ply"),
# scallops[,,1], ref, color=NULL, centered=T)
#PCA.scallop <- gm.prcomp(scallops)
#PC.sc <- PCA.scallop$x[,1]
#sc.preds <- shape.predictor(scallops, x= PC.sc, Intercept = FALSE,
# pred1 = min(PC.sc), pred2 = max(PC.sc)) # PC 1 extremes, more technically
#plotRefToTarget(ref, sc.preds$pred1)
#plotRefToTarget(ref, sc.preds$pred2)
#plotRefToTarget(ref, sc.preds$pred1, mesh = refmesh, method = "surface", mag = 1)
#plotRefToTarget(ref, sc.preds$pred2, mesh = refmesh, method = "surface", mag = 1)