A phylogenetic lollipop plot for phytools

Last night I tweeted about a phylogenetic “pin” or “lollipop” plot, inspired by some really nice visuals in a recent article by Pascal Title et al. (2024) & published in the journal Science.

Inspired by a plot style I learned from a new paper by @pascal_title et al. -- trait values plotted at the tips of a tree using "pins" in #Rstats #phytools with base R graphics: https://t.co/w9Xo1wkV6Z. pic.twitter.com/bcP031Tbu1

— Liam Revell (@phytools_liam) February 23, 2024

The tweet has been so popular, that I thought it would be worth putting into a function to be added (at some point) to phytools.

Here is that function. It involves a bunch of different tricks in base R & R graphics that I won’t get into – but interest R folks can take a look.

plotTree.lollipop<-function(tree,x,
  args.plotTree=list(),args.lollipop=list(),...){
  if(!inherits(x,c("matrix","data.frame"))) x<-as.matrix(x)
  h<-max(nodeHeights(tree))
  if(hasArg(panel_height)) panel_height<-list(...)$panel_height
  else panel_height<-1.0
  panel_height<-panel_height*h
  args.plotTree$tree<-tree
  args.plotTree$direction<-"upwards"
  if(is.null(args.plotTree$mar)) 
    args.plotTree$mar<-c(0.1,5.1,0.1,0.1)
  if(is.null(args.plotTree$ylim)) 
    args.plotTree$ylim<-c(0,h+ncol(x)*panel_height)
  if(is.null(args.plotTree$ftype)) 
    args.plotTree$ftype<-"off"
  if(is.null(args.plotTree$lwd)) args.plotTree$lwd<-1
  do.call(plotTree,args.plotTree)
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  if(pp$font){
    dx<-abs(diff(pp$x.lim))
    pdin<-par()$din[2]
    sh<-(pp$cex*strwidth(paste(" ",tree$tip.label,sep=""))+
      0.3*pp$cex*strwidth("W"))*(par()$din[1]/par()$din[2])*
      (diff(par()$usr[3:4])/diff(par()$usr[1:2]))
    new_h<-h+max(sh)
    panel_height<-(h-new_h+ncol(x)*panel_height)/ncol(x)
    h<-new_h
  }
  if(hasArg(ylab)) ylab<-list(...)$ylab
  else ylab<-if(!is.null(colnames(x))) 
    colnames(x) else rep("",ncol(x))
  for(i in ncol(x):1){
    d<-max(c(diff(range(x[,i])),max(x[,i])))
    y<-setNames(x[,i]/d*0.8*panel_height,rownames(x))
    lower<-h+(i-1)*panel_height+panel_height*0.05
    upper<-h+(i-1)*panel_height+panel_height*0.95
    polygon(c(0,max(pp$xx)+1,max(pp$xx)+1,0),
      c(lower,lower,upper,upper),
      border=FALSE,col="#F2F2F2")
    hh<-lower-min(c(0,min(y)))+0.05*panel_height
    lines(range(pp$xx),rep(hh,2),col="black",lty="dotted")
    segments(x0=pp$xx[1:Ntip(tree)],y0=rep(hh,Ntip(tree)),
      x1=pp$xx[1:Ntip(tree)],y1=y[tree$tip.label]+hh)
    labs<-pretty(c(min(c(0,min(x[,i]))),x[,i]),n=4)
    labs[!(labs>max(x[,i]))]->labs
    labs[!(labs<min(c(0,min(x[,i]))))]->labs
    axis(2,at=hh+max(y)/max(x[,i])*labs,
      labels=labs,las=1,cex.axis=0.6)
    args.lollipop$bg<-setNames(
      viridisLite::viridis(n=100)[ceiling(99*((y-
          min(y))/diff(range(y))))+1],
      names(y))
    args.lollipop$bg<-args.lollipop$bg[tree$tip.label]
    if(is.null(args.lollipop$pch)) args.lollipop$pch<-21
    if(is.null(args.lollipop$cex)) args.lollipop$cex<-1.2
    args.lollipop$x<-pp$xx[1:Ntip(tree)]
    args.lollipop$y<-y[tree$tip.label]+hh
    do.call(points,args.lollipop)
    mtext(ylab[i],2,line=3,at=mean(hh+max(y)/
      max(x[,i])*labs),cex=0.8)
  }
}

Let’s load phytools and try it out. This also requires installation of the package viridisLite.

library(phytools)

We can start by reiterating the example from last night’s post, but I’ll add a few more characters and show the tip labels of the tree, as follows.

data(anoletree)
anole_tree<-as.phylo(anoletree)
data(anole.data)
anole_data<-cbind(phyl.resid(anole_tree,
  x=as.matrix(anole.data[,"SVL",drop=FALSE]),
  Y=as.matrix(anole.data[,c(6,4,2)]))$resid,
  exp(anole.data[,"SVL",drop=FALSE]))
head(anole_data)

##                      TL        FLL          HL      SVL
## ahli        -0.11277698 0.19249584 -0.01872421 56.77693
## allogus     -0.09594253 0.18468706 -0.04134903 56.83430
## rubribarbus -0.09681106 0.16034890 -0.04588292 59.05505
## imias       -0.02764799 0.19445349 -0.10810397 60.32159
## sagrei       0.12995898 0.05913699 -0.09384291 58.39090
## bremeri      0.03559510 0.07113015 -0.11406836 61.15245

plotTree.lollipop(anole_tree,anole_data,
  args.plotTree=list(fsize=0.5,ftype="i"))

Awesome. For fun, let’s try body and home range size in mammals – a classic dataset from Garland et al. (1992).

data(mammal.tree)
data(mammal.data)
plotTree.lollipop(mammal.tree,log(mammal.data),
  ylab=c("log(body mass)","log(home range)"),
  args.plotTree=list(fsize=0.7,ftype="i"),
  args.lollipop=list(cex=1.5),
  panel_height=2.5)

Cool. Here’s some simulated data for Twitter.

tree<-ladderize(pbtree(n=200))
X<-fastBM(tree,nsim=6)
plotTree.lollipop(tree,X,
  args.plotTree=list(mar=c(1.1,2.1,1.1,1.1)),
  args.lollipop=list(cex=0.8))