I recently fielded a question about how, having fit an M*k* 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 M*k* 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)
```