Today I received the following email:

*“I'm a postdoc working on a project to reconstruct the ancestral … with … at the
University of …. We've been using phytools and getting nice results but we've run into
one issue we were hoping you might be able to help us with. We're interesting in generating
a posterior density map for our data, but, as far as we can tell, the function
densityMap can only be used for binary characters. Do you know of a way to do
a densityplot with non-binary characters (3 character states) or can you recommend a
workaround? Any assistance you can provide would be greatly appreciated. Thank you
for your time and your help.”*

This sounds like something I would *love* to implement, so perhaps I will; however
in the meantime there is a kind of workaround that can be done using `densityTree`

(not `densityMap`

, in fact) in the phytools package as follows:

```
library(phytools)
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 a b b b b a a c c a c a b a c c c c b c c c a a a
## Levels: a b c
```

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

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

```
## Done.
```

```
trees
```

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

```
colors<-setNames(c(rgb(1,0,0),rgb(0,1,0),rgb(0,0,1)),
letters[1:3])
densityTree(trees,method="plotSimmap",lwd=8,nodes="intermediate",
colors=colors,ylim=c(-5,26),compute.consensus=FALSE)
```

OK, so now we're part of the way there - but what about our legend?

Well, the following is a 'hack' - and would have to be adapted in its coordinates & so on for any particular plot - but hopefully you can see the idea:

```
densityTree(trees,method="plotSimmap",lwd=8,nodes="intermediate",
colors=colors,ylim=c(-5,26),compute.consensus=FALSE)
## make a triangle of cell coordinates
cells.x<-cells.y<-vector()
k<-1
asp<-abs(diff(par()$usr[4:3])/diff(par()$usr[2:1]))*par()$pin[2]/par()$pin[1]
for(i in 1:50){
for(j in i:(100-i+1)){
cells.x[k]<-0.65+(j-1)*0.3/100
cells.y[k]<-(-5)+(i-1)*2*0.3*asp/100
k<-k+1
}
}
## draw points
dx<-diff(range(cells.x))
dy<-diff(range(cells.y))
red.x<-max(cells.x)
red.y<-min(cells.y)
red<-1-sqrt((((cells.x-red.x)/dx)^2+((cells.y-red.y)/dy)^2)/2)
green.x<-min(cells.x)
green.y<-min(cells.y)
green<-1-sqrt((((cells.x-green.x)/dx)^2+((cells.y-green.y)/dy)^2)/2)
blue.x<-min(cells.x)
blue.y<-max(cells.y)
blue<-1-sqrt((((cells.x-blue.x)/dx)^2+((cells.y-blue.y)/dy)^2)/2)
points(cells.x,cells.y,col=rgb(red/max(red),green/max(green),blue/max(blue)),
pch=15,cex=0.5)
text(0.95,-5,"a",pos=2)
text(0.65,-5,"b",pos=4)
text((0.95-0.65)/2+0.65,max(cells.y),"c",pos=3)
```

OK, it's not perfect - but perhaps something to work from!

More later.

The tree & data were simulated as follows:

```
tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
Q<-matrix(c(-1,0.5,0.5,0.5,-1,0.5,0.5,0.5,-1),
3,3,dimnames=list(letters[1:3],letters[1:3]))
x<-as.factor(sim.history(tree,Q)$states)
```

## No comments:

## Post a Comment