1. Converting "phyl.pca" object to "princomp"

Various readers of the book have expressed interest in passing objects of class "phyl.pca" from the phytools function phyl.pca (Revell 2009) to visualization methods of other packages, such as factoextra. This is done easily using the phytools method as.princomp. Here’s a simple demo that follows the example of the book.

First, we can start by loading the R libraries that we’ll need. These include phytools (and its dependencies) and factoextra (which we will need to install from CRAN.

library(phytools)
library(factoextra)

Next, let’s load the tree and data for our analysis. In the book we used data files from this website, but these data are also packaged with phytools, so it will be easier for us to get them there.

data(anoletree)
data(anole.data)

To use factoextra in an interesting way we’re also going to need a grouping variable. Once again, in the book this was a separate input file; however, we happen to have our grouping variable mapped onto our tree in anoletree, so let’s pull them out using the phytools function getStates.

ecomorph<-getStates(anoletree,"tips")

This is what that grouping variable looks like. It encodes microhabitat specialization… e.g., TG = “trunk-ground”; GB = “grass-bush”; etc.

head(ecomorph,10)
##            ahli         allogus     rubribarbus           imias          sagrei         bremeri 
##            "TG"            "TG"            "TG"            "TG"            "TG"            "TG" 
## quadriocellifer      ophiolepis         mestrei           jubar 
##            "TG"            "GB"            "TG"            "TG"

Next, we’re ready for our phylogenetic principal components

ecomorph.pca<-phyl.pca(anoletree,anole.data)
ecomorph.pca
## Phylogenetic pca
## Standard deviations:
##        PC1        PC2        PC3        PC4        PC5        PC6 
## 0.33222579 0.09207406 0.05012082 0.04318123 0.02008655 0.01507125 
## Loads:
##            PC1         PC2         PC3         PC4         PC5         PC6
## SVL -0.9712234  0.16067288  0.01972570  0.14782215 -0.06211906  0.06935433
## HL  -0.9645111  0.16955087 -0.01203113  0.17994634  0.08064324 -0.04406887
## HLL -0.9814164 -0.02674808  0.10315533 -0.13790763  0.06887922  0.04126248
## FLL -0.9712265  0.17585694  0.10697935 -0.09105747 -0.06075142 -0.04864769
## LAM -0.7810052  0.37429334 -0.47398703 -0.15871456  0.00217418  0.00875408
## TL  -0.9014509 -0.42528918 -0.07614571  0.01709649 -0.01750404 -0.01088743

The object that comes out is of class "phyl.pca" – a special object class created for phytools. To use factoextra we just need to convert it to be a "princomp" object. phytools makes this easy using the as.princomp method as follows.

ecomorph.princomp<-as.princomp(ecomorph.pca)
ecomorph.princomp
## Call:
## ml_phyl.pca(tree = tree, Y = Y, method = method, mode = mode)
## 
## Standard deviations:
##     Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6 
## 0.33222579 0.09207406 0.05012082 0.04318123 0.02008655 0.01507125 
## 
##  6  variables and  observations.

Terrific. First, let’s try a screeplot using factoextra::fviz_screeplot as follows.

fviz_screeplot(ecomorph.princomp,addlabels=TRUE)
Screeplot from phylogenetic PCA of section 1.10 of the book.

Screeplot from phylogenetic PCA of section 1.10 of the book.

Lastly, let’s visualize the scores from our phylogenetic PCA, but with centroid plotted over each of our groups. To do that, we’ll use fviz_pca_ind as follows.

fviz_pca_ind(ecomorph.princomp,label="none",habillage=ecomorph,
  addEllipses=TRUE)
Scores for PC 1 and 2 of a phylogenetic PCA of anole morphology, with centroids and ellipses indicating different ecomorphological groupings of tip taxa.

Scores for PC 1 and 2 of a phylogenetic PCA of anole morphology, with centroids and ellipses indicating different ecomorphological groupings of tip taxa.

Neat!