Monday, April 24, 2017

Fixed, updated, & extended plotting method for "changesMap" object from density.multiSimmap method

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)

plot of chunk unnamed-chunk-2

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)

plot of chunk unnamed-chunk-2

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)

plot of chunk unnamed-chunk-3

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