Yesterday
I described an updated S3 density method for phytools' "multiSimmap"
object class.
Unfortunately, this update breaks the plotting method that also exists for the object
class (a "changesMap" object) typically returned by this method in its default
mode. This is because the plot method, though resulting in a nice visual, only
permitted a two-state (binary) character with two types of changes between states.
The following method fixes this issue by using different methods depending on whether the mapped character is binary or multistate:
plot.changesMap<-function(x,...){
p<-x$p
hpd<-x$hpd
bw<-x$bw
if(length(x$trans)==2){
plot(p[[1]]$mids,p[[1]]$density,xlim=c(min(x$mins)-1,
max(x$maxs)+1),ylim=c(0,1.2*max(c(p[[1]]$density,
p[[2]]$density))),
type="n",xlab="number of changes",
ylab="relative frequency across stochastic maps")
y2<-rep(p[[1]]$density,each=2)
y2<-y2[-length(y2)]
x2<-rep(p[[1]]$mids-bw/2,each=2)[-1]
x3<-c(min(x2),x2,max(x2))
y3<-c(0,y2,0)
polygon(x3,y3,col=make.transparent("red",0.3),border=FALSE)
lines(p[[1]]$mids-bw/2,p[[1]]$density,type="s")
y2<-rep(p[[2]]$density,each=2)
y2<-y2[-length(y2)]
x2<-rep(p[[2]]$mids-bw/2,each=2)[-1]
x3<-c(min(x2),x2,max(x2))
y3<-c(0,y2,0)
polygon(x3,y3,col=make.transparent("blue",0.3),border=FALSE)
lines(p[[2]]$mids-bw/2,p[[2]]$density,type="s")
add.simmap.legend(colors=setNames(c(make.transparent("red",0.3),
make.transparent("blue",0.3)),x$trans[1:2]),
prompt=FALSE,x=min(x$mins),y=0.95*par()$usr[4])
dd<-0.01*diff(par()$usr[3:4])
lines(hpd[[1]],rep(max(p[[1]]$density)+dd,2))
lines(rep(hpd[[1]][1],2),c(max(p[[1]]$density)+dd,
max(p[[1]]$density)+dd-0.005))
lines(rep(hpd[[1]][2],2),c(max(p[[1]]$density)+dd,
max(p[[1]]$density)+dd-0.005))
text(mean(hpd[[1]]),max(p[[1]]$density)+dd,
paste("HPD(",x$trans[1],")",sep=""),
pos=3)
lines(hpd[[2]],rep(max(p[[2]]$density)+dd,2))
lines(rep(hpd[[2]][1],2),c(max(p[[2]]$density)+dd,
max(p[[2]]$density)+dd-0.005))
lines(rep(hpd[[2]][2],2),c(max(p[[2]]$density)+dd,
max(p[[2]]$density)+dd-0.005))
text(mean(hpd[[2]]),max(p[[2]]$density)+dd,
paste("HPD(",x$trans[2],")",sep=""),
pos=3)
} else {
k<-length(x$states)
par(mfrow=c(k,k))
ii<-1
max.d<-max(unlist(lapply(p,function(x) x$density)))
for(i in 1:k){
for(j in 1:k){
if(i==j) plot.new()
else {
plot(p[[ii]]$mids,p[[ii]]$density,xlim=c(min(x$mins)-1,
max(x$maxs)+1),ylim=c(0,1.2*max.d),
type="n",xlab="number of changes",
ylab="relative frequency",main=x$trans[ii],font.main=1)
##title(main=)
y2<-rep(p[[ii]]$density,each=2)
y2<-y2[-length(y2)]
x2<-rep(p[[ii]]$mids-bw/2,each=2)[-1]
x3<-c(min(x2),x2,max(x2))
y3<-c(0,y2,0)
polygon(x3,y3,col=make.transparent("blue",0.3),border=FALSE)
lines(p[[ii]]$mids-bw/2,p[[ii]]$density,type="s")
dd<-0.03*diff(par()$usr[3:4])
lines(hpd[[ii]],rep(max(p[[ii]]$density)+dd,2))
text(mean(hpd[[ii]]),max(p[[ii]]$density)+dd,"HPD",pos=3)
ii<-ii+1
}
}
}
}
}
Here's a demo of what I mean:
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
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 1
## Levels: 0 1
y
## 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
## c a a a a a a a a c c c b b b b b c c c c c c a b a
## Levels: a b c
trees.x<-make.simmap(tree,x,nsim=100,model="ARD")
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## 0 1
## 0 -1.27244513 1.27244513
## 1 0.07744615 -0.07744615
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## 0 1
## 0.5 0.5
## Done.
trees.x
## 100 phylogenetic trees with mapped discrete characters
obj<-density(trees.x)
obj
##
## Distribution of changes from stochastic mapping:
## 0->1 1->0
## Min. :1 Min. :0
## Median :4 Median :1.5
## Mean :3.88 Mean :1.55
## Max. :8 Max. :4
##
## 95% HPD interval(0->1): [1, 7]
## 95% HPD interval(1->0): [0, 3]
plot(obj)
model<-matrix(
c(0,1,0,
1,0,1,
0,1,0),3,3)
rownames(model)<-colnames(model)<-letters[1:3]
model
## a b c
## a 0 1 0
## b 1 0 1
## c 0 1 0
trees.y<-make.simmap(tree,y,model=model,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## a b c
## a -0.3603485 0.3603485 0.0000000
## b 0.3603485 -0.7206969 0.3603485
## c 0.0000000 0.3603485 -0.3603485
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## a b c
## 0.3333333 0.3333333 0.3333333
## Done.
obj<-density(trees.y)
obj
##
## Distribution of changes from stochastic mapping:
##
## a->b a->c
## Min. :0 Min. :0
## Median :1 Median :0
## Mean :1.29 Mean :0
## Max. :5 Max. :0
##
## b->a b->c
## Min. :1 Min. :0
## Median :3 Median :4
## Mean :3.33 Mean :3.84
## Max. :7 Max. :6
##
## c->a c->b
## Min. :0 Min. :0
## Median :0 Median :1
## Mean :0 Mean :1.81
## Max. :0 Max. :6
##
## 95% HPD interval(a->b): [0, 4]
## 95% HPD interval(a->c): [0, 0]
## 95% HPD interval(b->a): [1, 5]
## 95% HPD interval(b->c): [2, 6]
## 95% HPD interval(c->a): [0, 0]
## 95% HPD interval(c->b): [0, 5]
plot(obj)
We can push the method to its extreme - although the figure may get a little messy:
data(anoletree)
ecomorph<-getStates(anoletree,"tips")
anoletrees<-make.simmap(anoletree,ecomorph,model="ER",nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## CG GB TC TG Tr Tw
## CG -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145
## GB 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145
## TC 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145
## TG 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145
## Tr 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145
## Tw 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## CG GB TC TG Tr Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
obj<-density(anoletrees)
obj
##
## Distribution of changes from stochastic mapping:
##
## CG->GB CG->TC
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.22 Mean :0.28
## Max. :2 Max. :2
##
## CG->TG CG->Tr
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.13 Mean :0.08
## Max. :1 Max. :1
##
## CG->Tw GB->CG
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.23 Mean :0.44
## Max. :3 Max. :3
##
## GB->TC GB->TG
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.63 Mean :0.95
## Max. :4 Max. :5
##
## GB->Tr GB->Tw
## Min. :0 Min. :0
## Median :0 Median :1
## Mean :0.25 Mean :1.08
## Max. :2 Max. :5
##
## TC->CG TC->GB
## Min. :0 Min. :0
## Median :1 Median :0
## Mean :1.27 Mean :0.71
## Max. :3 Max. :3
##
## TC->TG TC->Tr
## Min. :0 Min. :0
## Median :0 Median :1
## Mean :0.38 Mean :0.7
## Max. :3 Max. :3
##
## TC->Tw TG->CG
## Min. :0 Min. :0
## Median :1 Median :1
## Mean :1.05 Mean :0.99
## Max. :5 Max. :3
##
## TG->GB TG->TC
## Min. :0 Min. :0
## Median :3 Median :2
## Mean :3.24 Mean :1.92
## Max. :6 Max. :5
##
## TG->Tr TG->Tw
## Min. :0 Min. :0
## Median :1 Median :1
## Mean :0.96 Mean :1.59
## Max. :4 Max. :5
##
## Tr->CG Tr->GB
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.11 Mean :0.15
## Max. :2 Max. :3
##
## Tr->TC Tr->TG
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.18 Mean :0.17
## Max. :2 Max. :2
##
## Tr->Tw Tw->CG
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :0.21 Mean :0.71
## Max. :2 Max. :4
##
## Tw->GB Tw->TC
## Min. :0 Min. :0
## Median :1 Median :1
## Mean :1.5 Mean :1.6
## Max. :5 Max. :4
##
## Tw->TG Tw->Tr
## Min. :0 Min. :0
## Median :0 Median :0
## Mean :1.05 Mean :0.54
## Max. :5 Max. :3
##
## 95% HPD interval(CG->GB): [0, 1]
## 95% HPD interval(CG->TC): [0, 2]
## 95% HPD interval(CG->TG): [0, 1]
## 95% HPD interval(CG->Tr): [0, 1]
## 95% HPD interval(CG->Tw): [0, 1]
## 95% HPD interval(GB->CG): [0, 2]
## 95% HPD interval(GB->TC): [0, 3]
## 95% HPD interval(GB->TG): [0, 4]
## 95% HPD interval(GB->Tr): [0, 1]
## 95% HPD interval(GB->Tw): [0, 4]
## 95% HPD interval(TC->CG): [0, 3]
## 95% HPD interval(TC->GB): [0, 3]
## 95% HPD interval(TC->TG): [0, 2]
## 95% HPD interval(TC->Tr): [0, 2]
## 95% HPD interval(TC->Tw): [0, 4]
## 95% HPD interval(TG->CG): [0, 2]
## 95% HPD interval(TG->GB): [1, 6]
## 95% HPD interval(TG->TC): [0, 4]
## 95% HPD interval(TG->Tr): [0, 2]
## 95% HPD interval(TG->Tw): [0, 5]
## 95% HPD interval(Tr->CG): [0, 1]
## 95% HPD interval(Tr->GB): [0, 1]
## 95% HPD interval(Tr->TC): [0, 1]
## 95% HPD interval(Tr->TG): [0, 1]
## 95% HPD interval(Tr->Tw): [0, 1]
## 95% HPD interval(Tw->CG): [0, 3]
## 95% HPD interval(Tw->GB): [0, 4]
## 95% HPD interval(Tw->TC): [0, 4]
## 95% HPD interval(Tw->TG): [0, 4]
## 95% HPD interval(Tw->Tr): [0, 2]
plot(obj)
The tree & data for this exercise were simulated as follows:
tree<-pbtree(n=26,tip.label=LETTERS)
Q1<-matrix(c(-1,1,0.01,-0.01),2,2)
rownames(Q1)<-colnames(Q1)<-0:1
x<-as.factor(getStates(sim.history(tree,Q1),"tips"))
Q2<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3)
rownames(Q2)<-colnames(Q2)<-letters[1:3]
y<-as.factor(getStates(sim.history(tree,Q2),"tips"))
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.