I recently fielded a question about how, having fit an Mk discrete character evolution model and recontructed ancestral states, one might go about using this reconstruction to study variation in the rate of trait evolution from branch to branch in the phylogeny.
In point of fact, when we fit the Mk model in the first place, we had assumed that the evolutionary rate (or, specifically, our transition matrix Q) was homogeneous (that is, unchanging) across all the edges of the phylogeny.
Still, I believe that it is nonetheless quite reasonable to investigate the distribution of state changes across the tree – and one very handy technique for this is the method of stochastic character mapping (Huelsenbeck et al. 2003).
To illustrate what I mean by this, I will use two different datasets – one for eels (from Collar et al. 2014), and a second for primates (from Kirk & Kay 2004). Both can be obtained from the website of my recent Princeton University Press book with Luke Harmon.
Let’s start with the eel tree & dataset.
library(phytools)
eel.tree<-read.tree(file="http://www.phytools.org/Rbook/8/elopomorph.tre")
eel.data<-read.csv(file="http://www.phytools.org/Rbook/8/elopomorph.csv",
row.names=1,stringsAsFactors=TRUE)
head(eel.data)
## feed_mode Max_TL_cm
## Albula_vulpes suction 104
## Anguilla_anguilla suction 50
## Anguilla_bicolor suction 120
## Anguilla_japonica suction 150
## Anguilla_rostrata suction 152
## Ariosoma_anago suction 60
We’re going to analyze the trait feed_mode
: a binary ‘feeding mode’ character.
fmode<-setNames(eel.data$feed_mode,rownames(eel.data))
levels(fmode)
## [1] "bite" "suction"
Now let’s run stochastic mapping. In this case, from exercises in our
book,
I happen to recall that the "ER"
(equal-rates) model is the one that is
best supported by the data.
fmode.trees<-make.simmap(eel.tree,fmode,model="ER",nsim=1000,
pi="fitzjohn")
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## bite suction
## bite -0.01569437 0.01569437
## suction 0.01569437 -0.01569437
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5491439 0.4508561
## Done.
fmode.trees
## 1000 phylogenetic trees with mapped discrete characters
Normally, we might be concerned with the probability that each edge was in
each of the two mapped character states, and, if we were, we could easily
visualize that using the phytools method densityMap
. In this case we
will follow two other recent posts to this blog
(1,
2)
and re-color our plot using the ‘viridis’ color palette.
dMap.eels<-densityMap(fmode.trees,plot=FALSE,res=500)
## sorry - this might take a while; please be patient
pal<-viridisLite::viridis(n=20)
dMap.eels<-setMap(dMap.eels,pal)
plot(dMap.eels,fsize=c(0.6,0.9),lwd=4,legend=50)
obj<-summary(fmode.trees)
nn<-which(apply(obj$ace,1,max)<0.95)
nodelabels(node=as.numeric(names(nn)),
pie=obj$ace[nn,],piecol=pal[c(1,20)],cex=0.6)