PCA 3D

---

Contact: Email: thomas.puschelrouliez@anthro.ox.ac.uk Web: thomas.puschel.com

---

Yesterday I received an email asking me how did I generate the 3D plots displayed in this preprint so I decided to share the process in this post. For this example I used the data from:

Püschel T. A., Marcé-Nogué J., Gladman J. T., Bobe R., & Sellers W. I. (2018). Inferring locomotor behaviours in Miocene New World monkeys using finite element analysis, geometric morphometrics and machine-learning classification techniques applied to talar morphology. Journal of The Royal Society Interface, 15(146), 20180520.

#1) Load the required packages

library('pca3d')
library('geomorph')
# for the colour palette
library('NineteenEightyR')

#2) Load the data

load('Talar_shape.RData')

#4) PCA biplot

PCA_morpho = prcomp(z, scale=F)
eig_m<-(PCA_morpho$sdev)^2
eig1_m<-(PCA_morpho$sdev[1])^2
eig2_m<-(PCA_morpho$sdev[2])^2
variance1_m<-eig1_m*100/sum(eig_m)
variance2_m<-eig2_m*100/sum(eig_m)
xlab <- paste('Principal Component 1', '(', round(variance1_m,digits=2),  '%)')
ylab <- paste('Principal Component 2', '(', round(variance2_m,digits=2),  '%)')
# plot PC1 vs PC2
plot(PCA_morpho$x[,1], PCA_morpho$x[,2], pch=tmp$pch, cex=2.5, bg = tmp$colors, xlab=xlab, ylab=ylab, asp=T)#
legend(-0.1,-0.04, bty = "n",legend= unique(tmp$Locomotion), pch=unique(tmp$pch), col=unique(tmp$colors), pt.bg=unique(tmp$colors), cex=0.5)

#5) Interactive PCA 3D

Additional information about the rgl package can be found here

SH<-as.vector(c("sphere","cube","octahedron","tetrahedron"))
tmp$sh<-SH[as.numeric(as.factor(tmp$Locomotion))]

# Setup the background color
rgl.bg(color = "white") 

#pca 3d
pca3d(pca=PCA_morpho, components = 1:3, col = tmp$colors, title = NULL, 
      new = FALSE,axes.color = "black", bg = "white", radius = 3.5, 
      group = tmp$Locomotion,shape=tmp$sh, show.axes = TRUE,show.axe.titles = TRUE,
      axe.titles = NULL, show.plane = F,show.shadows = F, show.centroids = F,show.scale = T)

## [1] 0.008803568 0.004052094 0.004235430

# Groups
groups <- as.factor(tmp$Locomotion)
levs <- levels(groups)
group.col <- c(malibu(4))
#
x<-PCA_morpho$x[,1]
y<-PCA_morpho$x[,2]
z<-PCA_morpho$x[,3]

# Compute 95% CI ellipse for each one of the groups
for (i in 1:length(levs)) {
   group <- levs[i]
   selected <- groups == group
   xx <- x[selected]; yy <- y[selected]; zz <- z[selected]
   ellips <- ellipse3d(cov(cbind(xx,yy,zz)), 
                       centre=c(mean(xx), mean(yy), mean(zz)), level = 0.95) 
   shade3d(ellips, col = group.col[i], alpha = 0.1, lit = FALSE) 
   # show group labels
   # texts3d(mean(xx),mean(yy), mean(zz), text = group,
   #         col= group.col[i], cex = 2)
}
aspect3d(1,1,1)
rglwidget(webgl=TRUE, snapshot=FALSE)