General Introduction

In this lab session, and those to follow, we review the basics of how to conduct analyses in geomorph using RRPP. These tutorials may seem ‘thin’, but this is by design. We feel that R-coding is best learned by actually coding, and by a bit of directed trial and error. Thus our intention is to provide the direction, and you need to provide the initiative to make the analysis happen.

RECALL THAT ALL FUNCTIONS IN GEOMORPH AND RRPP HAVE RATHER EXTENSIVE HELP FILES. READ THEM!!!! They provide important details regarding the analysis, the input data, and most importantly, the various options of the function. There is considerable flexibility in one’s analysis when using geomorph and RRPP, so feel free to explore!

1: Reading Data

Before conducting any analysis one must first read one’s data into R. Base-R has many functions for reading datafiles of various types (e.g., read.csv, read.table, etc.) and geomorph adds to these by allowing morphometric-specific files to be read.

Important Functions

  • readland.nts
  • readland.tps
  • readland.shapes

Reading Data: Example

Here is a simple example for reading a tps style file:

mydata <- readland.tps("Data/salamanders.tps", specID ="imageID") # Specify specimen labels
## 
## No curves detected; all points appear to be fixed landmarks.
dim(mydata)
## [1]  12   2 287
mydata[,,1]
##            [,1]     [,2]
##  [1,]  8.818432 53.14680
##  [2,]  9.290848 52.18228
##  [3,]  5.235944 53.00901
##  [4,]  2.165240 51.47366
##  [5,]  1.909348 51.70987
##  [6,]  1.810928 51.80829
##  [7,]  1.102304 52.83186
##  [8,]  3.110072 54.48531
##  [9,]  6.023304 54.44594
## [10,]  8.916852 54.60342
## [11,] 10.649044 54.83962
## [12,] 14.015008 52.93028
NOTE: reading in landmarks and curves may be accomplished in several ways, including readland.shapes

Data Prep Components

At this step of the pipeline, one has some additional data ‘clean up’ to perform. For instance, if there are missing landamrks, one needs to estimate these (see Day 4 Materials). Likewise if there are positional differences due to articulations, or if one wishes to combine landmark configurations in a single analysis (e.g., heads and tails), one does that now (see Day 4 Materials).

2: Generalized Procrustes Analysis (GPA)

Here is a simple example of GPA superimposition

data(plethodon)
Y.gpa <- gpagen(plethodon$land, print.progress = F)

Plotting Specimens

View the data before and after GPA:

plotAllSpecimens(plethodon$land, links = plethodon$links)

plotAllSpecimens(Y.gpa$coords, links = plethodon$links)

Data Checking: Outliers

At this point, it is a good idea to check for outliers. If they are identified then one should go back to the original data and resolve any issues. Here is a simple example of how to check for these:

plotOutliers(Y.gpa$coords, inspect.outliers = T)

## 14 27  5 20  4 12 13 19 11 10  3  9  1 16  6 15  7 29  8 22  2 34 37 18 31 38 
## 14 27  5 20  4 12 13 19 11 10  3  9  1 16  6 15  7 29  8 22  2 34 37 18 31 38 
## 40 30 26 21 17 36 35 32 24 33 28 23 25 39 
## 40 30 26 21 17 36 35 32 24 33 28 23 25 39

3: Principal Components Analysis (PCA)

Here is a simple example of PCA

PCA <- gm.prcomp(Y.gpa$coords)
plot(PCA)

Maybe color some groups…

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 = unique(gps))

4: Thin-Plate Spline (TPS) Shape Visualization

Let’s visualize some shapes. First let’s look at a single specimen relative to the mean:

M <- mshape(Y.gpa$coords)

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))

Another option (if one has an outline)

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))
Advanced Shape Plotting: Shape Prediction

Predict a shape and plot this using shape.predictor (very useful!). This time along PC1:

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))

Visualizing 3D Shapes

One can deform a mesh along PC axes (or anything)

# (NOT RUN)
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)
plot(PCA.scallop, pch = 21, bg = "black", cex = 2)
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

plotRefToTarget(ref, sc.preds$pred1, mesh = refmesh, method = "surface", mag = 1)
plotRefToTarget(ref, sc.preds$pred2, mesh = refmesh, method = "surface", mag = 1)

# Saving deformed PLY
PC.example <- plotRefToTarget(ref, sc.preds$pred2, mesh = refmesh, method = "surface", mag = 1)
library(Rvcg)
vcgPlyWrite(PC.example, filename= "PC.example.ply", writeCol = FALSE)

5: Things to Explore on Your Own

  1. Read in your own data
  2. Check data for outliers (and fix challenging specimens)
  3. Perform GPA: visualize
  4. Perform PCA: visualize
  5. Obtain TPS visualizations
  6. Can you generate shapes for group means? For any specimen location in PC-space?
  7. Explore function options and output (there are many!)