Computing & visualizing the distribution of character transitions from a set of stochastic character mapped trees

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)