More on plotting multiple quantitative traits at the tips of the tree using bars

Yesterday, in response to a tweet, I posted a short demo illustrating how to simply create a plotTree.wBars fan-style graph in a semi-circle arc.

I think this came out pretty cool, so I thought I'd run another quick demo using a different dataset to illustrate the concepts.

Rather than provide a complete narrative description of what I'm doing, I have just commented my code a little bit below. Hopefully this will be enough for many readers to follow along.

The graph that results will illustrate the normalized (standardized to have the same variance) values of the first three PCs of a phylogenetic principal components analysis of a morphological dataset of anoles. The data come from Mahler et al. (2010) and also feature in my upcoming book with Luke Harmon.

## load packages
library(phytools)
## read data and tree from file
anole.data<-read.csv(file="http://www.phytools.org/Rbook/1/anole.data.csv",
    row.names=1)
anole.tree<-read.tree(file="http://www.phytools.org/Rbook/1/Anolis.tre")
## run phylogenetic PCA
anole.pca<-phyl.pca(anole.tree,anole.data,mode="corr")
anole.pca

## Phylogenetic pca
## Standard deviations:
##       PC1       PC2       PC3       PC4       PC5       PC6 
## 2.2813006 0.6678948 0.4673306 0.3180545 0.1408844 0.1008910 
## Loads:
##            PC1          PC2         PC3         PC4           PC5
## SVL -0.9763368 -0.015638858  0.18242357  0.07664922  0.0496637309
## HL  -0.9600799 -0.048076434  0.22750190  0.13558866 -0.0596373555
## HLL -0.9699480  0.137476262 -0.05120924 -0.17065492 -0.0871581061
## FLL -0.9770212 -0.001554872  0.07778124 -0.17827690  0.0753072554
## LAM -0.8233484 -0.511658499 -0.24300346  0.03536477 -0.0009539185
## TL  -0.8695250  0.403523205 -0.25620191  0.12142006  0.0235939736
##               PC6
## SVL  0.0700136069
## HL  -0.0472999428
## HLL  0.0309702558
## FLL -0.0438572618
## LAM  0.0008967491
## TL  -0.0125051491

## extract scores for first three PCs
anole.pc<-scores(anole.pca)[,1:3]
## normalize to have the same variance
anole.norm<-anole.pc/matrix(rep(apply(anole.pc,2,sd),nrow(anole.pc)),
    nrow(anole.pc),ncol(anole.pc),byrow=TRUE)*0.2
## set background black & foreground to white for plotting
par(fg="white",bg="black")
## compute max tree height and number of traits
h<-max(nodeHeights(anole.tree))
m<-ncol(anole.pc)
## set colors
cols<-RColorBrewer::brewer.pal(n=m,"Spectral")
## set x & y limits
xlim<-ylim<-1.2*c(-h,h)+c(-1,1)*0.15*m*h+0.2*c(-h,h)
ylim<-c(0,ylim[2])
## plot tree
plotTree(anole.tree,type="fan",xlim=xlim,ylim=ylim,
    lwd=1,ftype="i",fsize=0.5,part=0.5,color="white")
## add traits
for(i in 1:m){
    tt<-anole.tree ## copy tree
    ## extend tip edges
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
        tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+
        0.35*h+(i-1)*0.2*h
    ## plot transparently
    plotTree(tt,color="transparent",type="fan",xlim=xlim,
        ylim=ylim,lwd=1,ftype="off",add=TRUE,part=0.5)
    pp1<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    ## extend again
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
        tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+
        0.2*h
    ## clip plot area
    clip(x1=par()$usr[1],x2=par()$usr[2],y1=ylim[1],y2=par()$usr[4])
    ## draw circle
    plotrix::draw.circle(0,0,radius=h+0.35*h+(i-1)*0.2*h,
        border="#B8B8B8")
    ## unclip
    clip(x1=par()$usr[1],x2=par()$usr[2],y1=par()$usr[3],
        y2=par()$usr[4])
    ## plot transparently again
    plotTree(tt,color="transparent",type="fan",xlim=xlim,
        ylim=ylim,lwd=1,ftype="off",add=TRUE,part=0.5)
    pp2<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    ## graph lines
    par(lend=1)
    for(j in 1:Ntip(anole.tree)){
        ii<-which(rownames(anole.pc)==anole.tree$tip.label[j])
        dx<-(pp2$xx[j]-pp1$xx[j])*anole.norm[ii,i]
        dy<-(pp2$yy[j]-pp1$yy[j])*anole.norm[ii,i]
        lines(pp1$xx[j]+c(0,dx),pp1$yy[j]+c(0,dy),lwd=8,
            col=cols[i])
    }
}
## create custom legend
xx<-rep(0.85*par()$usr[4],m)
yy<-0.95*par()$usr[4]-1.5*1:m*strheight("W")
text(xx,yy,colnames(anole.norm),pos=4,cex=0.8)
scale.bar<-sqrt(diag(anole.pca$Eval)[1:3])*0.2*h
xx2<-xx-scale.bar
segments(x0=xx,y0=yy,x1=xx2,y1=yy,lwd=8,col=cols)
text(x=xx2[1],y=0.95*par()$usr[4],"standard deviations",
    pos=4,cex=0.8)

That's it!