Plotting a quantitative trait on the tree using variable branch widths: Part III

The other day I tweeted about a recent post on this blog in which I showed how to vary line widths both within & between edges of the tree as a function of the observed or reconstructed value of a quantitative trait.

As promised - here's how to do it using ape, #phytools, & #rstats. Mapping a continuous character onto the tree using variable edge widths: https://t.co/yHoscQWtff. https://t.co/l9KHZ6PWCH pic.twitter.com/hU2sAE4CRI

— Liam Revell (@phytools_liam) July 16, 2020

I did this by: (1) breaking up the edges of the tree into many singleton nodes; (2) reconstructing ancestral character states for the continuous trait at each and every one of these nodes; and (3) plotting the tree using ape::plot.phylo which permits different line widths for each edge.

Although this came out pretty nicely, something that bothered me is that by using a simple call to segments to vary our edge widths, fatter plotted lines of the tree also end up getting a tiny bit longer. How annoying!

Today, I set out to solve this - and I have, though in the end it is by creating a new function & set of methods that are decidedly less flexible than my previous, inarguably hacky, solution. How is it less flexible? Well, for one, since I have written a new plotting method, it will only work for right- and left-facing phylograms, and not for other tree plotting styles supported by plot.phylo or phytools::plotTree.

Here is my function code:

## function to create "edge.widthMap" object
edge.widthMap<-function(tree,x,res=100,...){
    if(!inherits(tree,"phylo")) 
        stop("tree should be an object of class \"phylo\".")
    tree<-as.phylo(tree)
    #h<-max(nodeHeights(tree))
    #LL<-setNames(seq(0,h,length.out=res),1:res)
    #tree<-map.to.singleton(make.era.map(tree,LL))
    a<-fastAnc(tree,x)
    node.values<-c(x[tree$tip.label],a)
    edge.values<-apply(tree$edge,1,function(e,nv)
        mean(nv[e]),nv=node.values)
    edge.widths<-edge.values
    object<-list(tree=tree,edge.widths=edge.widths,
        node.values=node.values)
    class(object)<-"edge.widthMap"
    object
}
## print method
print.edge.widthMap<-function(x,...){
    cat("Object of class \"edge.widthMap\" containing:\n")
    cat(paste("(1) Phylogenetic tree with",Ntip(x$tree),
        "tips and",x$tree$Nnode,"internal nodes.\n"))
    cat("(2) Vector of edge widths corresponding to the values of a mapped\n")
    cat("    quantitative trait.\n\n")
}
## plot method
plot.edge.widthMap<-function(x,max.width=1,legend="trait value",...){
    if(hasArg(min.width)) min.width<-list(...)$min.width
    else min.width<-0
    h<-max(nodeHeights(x$tree))
    node.values<-x$node.values-min(x$node.values)
    node.values<-node.values*((max.width-min.width)/
        max(node.values))+min.width
    args.list<-list(...)
    args.list$tree<-x$tree
    args.list$type<-"phylogram"
    if(!is.null(args.list$direction)){
        if(!args.list$direction%in%c("leftwards","rightwards"))
            args.list$direction<-"rightwards"
    } else args.list$direction<-"rightwards"
    if(is.null(args.list$ylim)) 
        args.list$ylim<-c(1,Ntip(x$tree)+Ntip(x$tree)/25)
    if(is.null(args.list$ftype)) args.list$ftype<-"i"
    if(is.null(args.list$fsize)) 
        args.list$fsize<-36*par()$pin[2]/par()$pin[1]/
            Ntip(x$tree)
    if(is.null(args.list$color)){ 
        args.list$color<-"transparent"
        color<-"gray62"
    } else {
        color<-args.list$color
        args.list$color<-"transparent"
    }
    do.call(plotTree,args.list)
    obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    asp<-par()$pin[2]/par()$pin[1]/
        (diff(par()$usr[4:3])/diff(par()$usr[2:1]))
    for(i in 1:nrow(x$tree$edge)){
        polygon(x=c(obj$xx[c(x$tree$edge[i,1],
            x$tree$edge[i,1:2],
            x$tree$edge[i,2:1],
            x$tree$edge[i,1])])+c(0,
            asp*node.values[x$tree$edge[i,1]]/2,
            0,0,asp*node.values[x$tree$edge[i,1]]/2,0),
            y=rep(obj$yy[x$tree$edge[i,2]],6)+
            c(node.values[x$tree$edge[i,1]],
            node.values[x$tree$edge[i,1:2]],
            -node.values[x$tree$edge[i,2:1]],
            -node.values[x$tree$edge[i,1]])/2,
            border=FALSE,col=color)
    }
    for(i in 1:x$tree$Nnode+Ntip(x$tree)){
        nn<-x$tree$edge[which(x$tree$edge[,1]==i),2]
        yy<-range(obj$yy[nn])
        polygon(x=rep(obj$xx[i],4)+
            asp*c(-node.values[i]/2,node.values[i]/2,
            node.values[i]/2,-node.values[i]/2),
            y=c(rep(yy[1],2),rep(yy[2],2))+
            c(-rep(node.values[i],2),
            rep(node.values[i],2))/2,
            border=FALSE,col=color)
    }
    leg.length<-0.4*h
    polygon(x=c(0,0,leg.length,leg.length),
        y=Ntip(x$tree)+Ntip(x$tree)/25+
        c(-min.width/2,min.width/2,max(node.values)/2,
        -max(node.values)/2),
        border=FALSE,col=color)
    text(0,Ntip(x$tree)+Ntip(x$tree)/25-0.1*max.width,
        round(min(x$node.values),2),pos=1,
        cex=0.8)
    text(leg.length,Ntip(x$tree)+Ntip(x$tree)/25-0.1*max.width,
        round(max(x$node.values),2),pos=1,cex=0.8)
    text(leg.length/2,Ntip(x$tree)+Ntip(x$tree)/25-0.1*max.width,
        legend,pos=1,cex=0.8)
}

Now, let's try it:

library(phytools)
data(sunfish.tree)
data(sunfish.data)
gape.width<-setNames(sunfish.data$gape.width,
    rownames(sunfish.data))
## create plotting object
gw.object<-edge.widthMap(sunfish.tree,gape.width)
plot(gw.object,fsize=0.9,max.width=0.8,
    legend="relative gape width")

data(anoletree)
data(anole.data)
svl<-setNames(exp(anole.data$SVL),rownames(anole.data))
svl.object<-edge.widthMap(anoletree,svl)
plot(svl.object,max.width=0.9,color=palette()[2],
    legend="SVL (mm)",min.width=0.05)

Awesome, right?

The hardest thing was figuring out how to get the corners right. I won't bore the reader by explaining how I did it - but if you're interested, you can look at the code & see!

Lastly, for this kind of plot with slanted lines - or shapes, it is absolutely key that you do not plot to the default graphical device in R, except for testing. I recommend using pdf( ) or png( ), the latter with res=600 or above. For example:

 png(file="anole-tree.png",width=7,height=9,units="in",
    res=1200)
plot(svl.object,max.width=0.9,color=palette()[2],
    legend="SVL (mm)",min.width=0.05)
dev.off()

That's it for now.