Generating a "contMap" visualization with node values from an early-burst model

I recently posted on reconstructing ancestral states for a continuously valued character under the 'early burst' (EB) trait evolution model.

However, the original user question also involved projecting these reconstructed values onto the tree visually using the phytools function, contMap. This can be done - kind off. What we must do, in fact, is obtain reconstructed node values using anc.ML, as I showed in my previous post, and then force these values to become the node states in contMap. The result will show color values at nodes reconstructed under EB - but values along edges interpolated assuming the traditional BM model. There are ways we could imagine getting around this, but they are much more complicated than what I'll demonstrate here.

First, let's start with our tree & data. This is the same tree & data that I used in my prior post:

library(phytools)
packageVersion("phytools")

## [1] '0.6.21'

tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  t29, t43, t44, t68, t75, t86, ...
## 
## Rooted; includes branch lengths.

x

##          t29          t43          t44          t68          t75 
##  0.438496287  0.472141210  0.508929413  0.476404420  0.500218520 
##          t86          t87          t56          t57          t21 
##  0.442388038  0.412293060  0.325040293  0.375942650  0.420729607 
##          t32          t33          t71          t72          t18 
##  0.144507151  0.074422614  0.188952192  0.203205398  0.156567264 
##           t9          t30          t48          t49          t37 
##  0.228356432  0.201488956  0.150009612  0.084055077  0.151670548 
##          t23          t24           t7           t5          t88 
##  0.120957147  0.276016226  0.289390081  0.286052945  0.285526927 
##          t89           t1          t58          t59          t84 
##  0.269553017  0.338785228 -0.638255514 -0.631750846 -0.655941964 
##          t85          t16          t65          t66          t95 
## -0.640326625 -0.817037017 -0.685259524 -0.663680995 -0.659381902 
##          t96          t90          t79          t22          t17 
## -0.676857787 -0.669947136 -0.673772868 -0.617029153 -0.709658998 
##          t82          t83           t8          t10          t19 
## -0.696291418 -0.688863550 -0.978188152 -0.876536606 -0.744679735 
##          t20           t4          t50          t51          t67 
## -0.782603431 -0.766160613 -0.779978124 -0.763065919 -0.877586463 
##          t93          t94          t91          t92          t27 
## -0.839376372 -0.853901636 -0.636389151 -0.643581270 -0.652135509 
##          t28           t6          t42          t80          t81 
## -0.657921272 -0.964749242 -0.994739494 -0.971873925 -0.970999109 
##          t41          t52          t53          t62          t78 
## -1.047142667 -1.030997961 -1.042163333 -1.037583807 -1.022104291 
##          t99         t100          t34          t35          t45 
## -1.025954239 -1.012997986 -1.047053377 -0.987320158 -1.086046548 
##          t46          t36          t38          t39          t54 
## -1.156241374 -1.120119673 -1.075409445 -1.036747839 -0.663847989 
##          t55          t97          t98          t47          t31 
## -0.698365634 -0.405954828 -0.414290095 -0.383287275 -0.344969991 
##          t11          t13          t76          t77          t73 
## -0.370850161 -1.095050380 -1.037155643 -1.045283871 -0.952923639 
##          t74           t2           t3          t14          t63 
## -1.001848791  0.141608367  0.061594482  0.313000952  0.091721162 
##          t64          t69          t70          t40          t60 
##  0.102012849  0.068840893  0.072838839  0.060826002  0.001699522 
##          t61          t26          t25          t15          t12 
##  0.065996008  0.018760874  0.174547937  0.060346979  0.097391382

Now, let's get our ML ancestral states under the EB model as follows:

fitEB<-anc.ML(tree,x,model="EB")
fitEB

## Ancestral character estimates using anc.ML under a(n) EB model:
##       101       102       103       104       105       106       107 
## -0.149645  0.303437  0.456207  0.462999  0.483923  0.463864  0.461351 
##       108       109       110       111       112       113       114 
##  0.441226  0.269318  0.254532  0.256832  0.341497  0.349491  0.134272 
##       115       116       117       118       119       120       121 
##  0.243730  0.188746  0.195718  0.237630  0.205890  0.170694  0.146584 
##       122       123       124       125       126       127       128 
##  0.128694  0.246328  0.247085  0.241900  0.216083  0.277122 -0.169551 
##       129       130       131       132       133       134       135 
## -0.682232 -0.681542 -0.688835 -0.691968 -0.647491 -0.639588 -0.648064 
##       136       137       138       139       140       141       142 
## -0.671030 -0.674188 -0.669582 -0.670976 -0.669138 -0.668707 -0.680070 
##       143       144       145       146       147       148       149 
## -0.786477 -0.812172 -0.799949 -0.794589 -0.800530 -0.764398 -0.694839 
##       150       151       152       153       154       155       156 
## -0.807828 -0.807780 -0.792406 -0.773237 -0.812107 -0.785678 -0.851186 
##       157       158       159       160       161       162       163 
## -0.847010 -0.688502 -0.640400 -0.664724 -0.917186 -0.977880 -0.972587 
##       164       165       166       167       168       169       170 
## -0.811511 -1.049309 -1.042513 -1.039931 -1.054948 -1.054235 -1.032098 
##       171       172       173       174       175       176       177 
## -1.025175 -1.019691 -1.054352 -1.043612 -1.100283 -1.067161 -1.064348 
##       178       179       180       181       182       183       184 
## -0.683323 -0.392338 -0.378467 -0.391821 -0.409841 -1.036606 -1.013423 
##       185       186       187       188       189       190       191 
## -1.037257 -0.983257  0.044576  0.095546  0.121583  0.166058  0.101472 
##       192       193       194       195       196       197       198 
##  0.106063  0.072185  0.094356  0.088166  0.084205  0.070487  0.052460 
##       199 
##  0.039961 
## 
## Fitted model parameters & likelihood:
##   sigma^2         r   logLik
##  0.449798 -4.169326 423.7552
## 
## R thinks it has found the ML solution.

Finally, let's generate our contMap visualization forcing the internal nodes to have color values determined by our EB reconstruction:

EB<-contMap(tree,x,method="user",anc.states=fitEB$ace,plot=FALSE)
EB

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
## 
## (2) A mapped continuous trait on the range (-1.156241, 0.508929).

plot(EB,ftype="off",lwd=c(3,6))

For fun, we can could compare this to a simpler BM model to see if we can detect a difference in the visualization:

BM<-contMap(tree,x,plot=FALSE)
BM

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
## 
## (2) A mapped continuous trait on the range (-1.156241, 0.508929).

plot(BM,ftype="off",lwd=c(3,6))

I'm color-blind, so you can be the judge.

Note that, as always, we have lots of user control over our contMap visualizations, for instance we can very easily change the color scheme of our plots as follows:

par(mfrow=c(1,2))
EB<-setMap(EB,topo.colors(n=10))
plot(EB,ftype="off",lwd=c(4,6),xlim=c(-0.1,1.1),legend=1,
    mar=c(0.1,0.1,3.1,0.1),outline=FALSE)
title(main="contMap under \"EB\" model")
BM<-setMap(BM,topo.colors(n=10))
plot(BM,ftype="off",lwd=c(4,6),legend=FALSE,xlim=c(-0.1,1.1),
    ylim=get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim,
    mar=c(0.1,0.1,3.1,0.1),outline=FALSE)
title(main="contMap under \"BM\" model")

It's a weird color scheme; however, (to my eye) it does help to visualize the larger change concentrated towards the root. Note that the data are the same in both case, so no matter our reconstruction we will tend to see the EB pattern.

Well, that's basically the idea.