Even more about plotting a discrete & continuous traits at the tips of a "fan" or "arc" style tree

Earlier today I tweeted about graphing discrete & continuous traits at the tips of a plotted tree in “arc” or “fan” style. The way I was doing this was just by added thickened lines that “grow” out of the tips of the tree.

A response from Ted Stankowich, however, has persuaded me to re-do this using polygons – even though the geometry is substantially more complicated. (It still requires nothing more than high school trig – but nonetheless!)

We've wrestled with whether to put the species names at the branch tips and the color rings outside of them. This plot is cool looking, but it seems like it's very hard to trace back the species name to its tip through the traits. Having white borders between markers might help? pic.twitter.com/uXdJNTh1wR

— Ted Stankowich 🦨🐆🦌🦔🦏🦡 (@CSULBMammalLab) March 27, 2024

The main advantages of polygons is that it allows us to incorporate “spacers” between adjacent species values, and easily widen the polygons as we go from the inside to the outside of our plotted traits.

OK. Here’s the code.

plotFanTree.wTraits<-function(tree,X,type=c("arc","fan"),
  ...){
  X<-if(is.vector(X)) as.matrix(X[tree$tip.label]) else 
    X[tree$tip.label,]
  h<-max(nodeHeights(tree))
  d<-min(ncol(X)*0.07*h,h)
  type<-type[1]
  if(!(type%in%c("arc","fan"))) type<-"fan"
  ftype<-if(hasArg(ftype)) list(...)$ftype else "i"
  fsize<-if(hasArg(fsize)) list(...)$fsize else 0.5
  part<-if(hasArg(part)) list(...)$part else 
    min(0.99,(Ntip(tree)-2)/Ntip(tree))
  arc_height<-if(hasArg(arc_height)) list(...)$arc_height else 
    0.7
  spacer<-if(hasArg(spacer)) list(...)$spacer else 0.025
  spacer<-spacer*(2*pi*part/(Ntip(tree)-1))/2
  xlim<-if(hasArg(xlim)) list(...)$xlim else NULL
  ylim<-if(hasArg(ylim)) list(...)$ylim else NULL
  if(hasArg(colors)) colors<-list(...)$colors
  else {
    colors<-list()
    for(i in 1:ncol(X)){
      if(is.numeric(X[,i])){
        colors[[i]]<-setNames(hcl.colors(n=100),1:100)
      } else {
        if(!is.factor(X[,i])) X[,i]<-as.factor(X[,i])
        colors[[i]]<-setNames(
          palette.colors(n=length(levels(X[,i]))),
          levels(X[,i]))
      }
    }
  }
  tt<-tree
  tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+d
  plotTree(tt,type=type,ftype=ftype,fsize=fsize,
    part=part,color="transparent",
    arc_height=arc_height*h/max(nodeHeights(tt)),
    xlim=xlim,ylim=ylim)
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  outer_rad<-max(pp$xx)
  plotTree(tree,type=type,ftype="off",part=part,
    lwd=1,add=TRUE,xlim=pp$x.lim,ylim=pp$y.lim,
    arc_height=arc_height,ftype="off")
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  inner_rad<-max(pp$xx)
  par(lend=3)
  for(i in 1:ncol(X)){
    if(is.numeric(X[,i])){
      x_seq<-seq(min(X[,i]),max(X[,i]),length.out=100)
      x_ind<-sapply(X[,i],function(x,y) which.min((x-y)^2),
        y=x_seq)
      colors[[i]]<-colorRampPalette(colors[[i]])(n=100)
      cols<-colors[[i]][x_ind]
    } else {
      cols<-colors[[i]][X[tree$tip.label,i]]  
    }
    for(j in 1:Ntip(tree)){
      start<-if(pp$xx[j]>0) 
        (i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)) else 
          -((i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)))
      end<-if(pp$xx[j]>0) i*d/ncol(X) else -i*d/ncol(X)
      th<-atan(pp$yy[j]/pp$xx[j])
      theta<-(2*pi*part/(Ntip(tree)-1))/2-spacer
      sign<-if(pp$xx[j]>0) 1 else -1
      H1<-(sign*inner_rad+start)/cos(theta)
      H2<-(sign*inner_rad+end)/cos(theta)
      th_up<-th+theta
      th_down<-th-theta
      x<-c(H1*cos(th_down),H2*cos(th_down),
        H2*cos(th_up),H1*cos(th_up))
      y<-c(H1*sin(th_down),H2*sin(th_down),
        H2*sin(th_up),H1*sin(th_up))
      polygon(x,y,col=cols[j],border=FALSE)
    }
  }
  invisible(colors)
}

First, let’s apply it to the example from earlier today to see that it still works! Here, I’ll set the optional argument spacer = 0 so that we have no white space between species levels of each character.

library(phytools)
data("liolaemid.data")
data("liolaemid.tree")
liolaemid.data[liolaemid.tree$tip.label,]->liolaemid.data
colors<-list(
  c("blue","white","red"),
  terrain.colors(n=10),
  setNames(c("#F0EAD6","#DF536B"),c("O","V")))
cols<-plotFanTree.wTraits(liolaemid.tree,
  liolaemid.data[,3:1],lwd=12,colors=colors,ftype="off",
  spacer=0)
legend(x=0,y=0.7*max(nodeHeights(liolaemid.tree)),
  names(colors[[3]]),lwd=8,col=colors[[3]],
  title="parity mode",bty="n",xjust=0.5,yjust=0.5)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[2]],
  title="maximum altitude (m)",
  lims=range(liolaemid.data[,2]),digits=2,prompt=FALSE,
  x=-0.75*max(nodeHeights(liolaemid.tree)),
  y=0.2*max(nodeHeights(liolaemid.tree)),subtitle="",
  lwd=8,outline=FALSE)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[1]],
  title="environmental temp.",
  lims=range(liolaemid.data[,3]),digits=2,prompt=FALSE,
  x=-0.75*max(nodeHeights(liolaemid.tree)),
  y=-0.15*max(nodeHeights(liolaemid.tree)),subtitle="",
  lwd=8,outline=FALSE)

So far, so good.

Now, let’s try a larger tree & more traits. Here, I’ll use a phylogeny of Mycalesina butterflies from Halali et al. (2020), but combined with 20 simulated continuous traits.

data("butterfly.tree")
X<-fastBM(butterfly.tree,nsim=20)
plotFanTree.wTraits(butterfly.tree,X,fsize=0.4)

Awesome. Let’s “zoom in” (actually, just re-plot but with xlim and ylim changed) so that we can see that we can see how the spacers look close up.

plotFanTree.wTraits(butterfly.tree,X,fsize=1,xlim=c(-80,-15),
  ylim=c(-50,-15),spacer=0.05)

That’s pretty much the idea, right?