## Saturday, July 18, 2020

### 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.

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)
``````