I had a query today about how to overlay the posterior probabilities from stochastic mapping on a single (“representative,” whatever that means) stochastic map tree. An example of this can be seen in Figure 4.3 of Revell (2014): my chapter on plotting methods in a recent edited volume on phylogenetic comparative methods.

It is not too hard. So, having loaded phytools:

```
library(phytools)
```

imagine we have a phylogeny (`tree`

) and a discrete character
vector `x`

:

```
tree
```

```
##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
##
## Rooted; includes branch lengths.
```

```
x
```

```
## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
## a b a a a a a b a a a a a a a b a a a b b a b b b b
## Levels: a b
```

Now, let's do stochastic mapping. We have multiple options for this. Here
I will use `model="ER"`

(equal rates) and `Q="empirical"`

(use the MLE of `Q`

):

```
mapped.trees<-make.simmap(tree,x,nsim=100,model="ER",Q="empirical")
```

```
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## a b
## a -0.008508737 0.008508737
## b 0.008508737 -0.008508737
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## a b
## 0.5 0.5
```

```
## Done.
```

```
mapped.trees
```

```
## 100 phylogenetic trees with mapped discrete characters
```

Now, let's compute the posterior probabilities at all the nodes of the tree.
We can do this using the S3 method `summary`

for objects of
class `"multiSimmap"`

:

```
obj<-summary(mapped.trees)
obj
```

```
## 100 trees with a mapped discrete character with states:
## a, b
##
## trees have 8.82 changes between states on average
##
## changes are of the following types:
## a,b b,a
## x->y 5.9 2.92
##
## mean total time spent in each state is:
## a b total
## raw 709.2151934 336.9705730 1046.186
## prop 0.6779056 0.3220944 1.000
```

The posterior probabilities at nodes are in `obj$ace`

and we
can plot these as follows (once again, using a custom S3 method):

```
colors<-setNames(c("blue","red"),c("a","b"))
plot(obj,colors=colors,ftype="i")
add.simmap.legend(colors=colors,prompt=FALSE,x=0,y=25)
```

Cool. We can also, however, overlay these posterior probabilities on a 'representative' stochastic map tree from the posterior distribution. Let's pick one such tree at random, plot it, and add the posterior probabilities to nodes:

```
i<-sample(1:length(mapped.trees),1)
plot(mapped.trees[[i]],colors=colors,lwd=4,ftype="i")
nodelabels(pie=obj$ace,piecol=colors[colnames(obj$ace)],cex=0.6)
```

That's it!

## No comments:

## Post a Comment