Reminder on how to overlay posterior probabilities from stochastic mapping on a single 'representative' stochastic map tree

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!

Phylogenetic regression when estimated Pagel's λ > 1

A phytools user reported to me the other day that they were finding that a PGLS model that was non-significant assuming a Brownian correlation structure of the residual error (corBrownian) became highly significant with Pagel's λ just slightly greater than 1.0.

This seems kind of mysterious because λ=1.0 is equivalent to an assumption of Brownian motion. On the other hand, λ is only defined when the correlations implied by λ are <=1.0. Usually this is only for values of λ (no more than) slightly larger than 1. [λ is always defined between 0 & 1.]

My hunch is that this can be best understood using contrasts. With contrasts, when we have a contrast between terminal species of zero or near zero length, the inevitable result for λ>1, this contrast (being standardized by dividing by a quantity very close to zero) will tend to have disproportionate weight in the contrasts regression.

To examine this, I set up the following simple simulation. First, I simulated a tree with slightly longer terminal edges than expected under pure-birth. This is probably not necessary, but I did this just to help exagerrate the effect of λ>1. Next, I simulated uncorrelated data on that tree assuming BM for x & y. Then, I fit a PGLS either setting λ to 1 (i.e., BM) or setting λ to it's maximum possible value. Finally I obtained the P-values for each fitted model & compared them.

Here is what I found:

## load libraries
library(phytools)
library(nlme)

## function for simulation
foo<-function(nrep=500){

    Pbm<-Pml<-lambda<-vector()

    cat("\n|")

    for(i in 1:nrep){
        tree<-phytools:::lambdaTree(pbtree(n=26,tip.label=LETTERS),0.95)

        x<-fastBM(tree)
        y<-fastBM(tree)

        corBM<-corBrownian(1,tree)
        lambda[i]<-phytools:::maxLambda(tree)
        corML<-corBrownian(1,phytools:::lambdaTree(tree,lambda[i]))

        fitbm<-gls(y~x,data=data.frame(x,y),correlation=corBM)
        fitml<-gls(y~x,data=data.frame(x,y),correlation=corML)

        Pbm[i]<-anova(fitbm)$'p-value'[2]
        Pml[i]<-anova(fitml)$'p-value'[2]

        cat(".")
        if(i%%50==0&&i!=nrep) cat("\n ")
        flush.console()

    }

    cat("|\n\n")

    cat("\nDone simulation.")

    list(Pbm=Pbm,Pml=Pml,lambda=lambda)

}
## run
obj<-foo()

## 
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
## Done simulation.

## visualize our results
par(mfrow=c(2,1))
hist(obj$Pbm,col="grey",xlab="P-value assuming BM",main="")
hist(obj$Pml,col="grey",xlab="P-value assuming lambda=max(lambda)",
    main="")

hist(obj$lambda,col="grey",main="Maximum value of lambda",xlab="lambda")

I'm not sure what to make of this - and whether this result suggests (or not) that λ should only be set or estimated between 0 & 1. This is already the strategy employed, I believe, by functions for PGLS in the caper package. Nonetheless, the result is here. Make of it what you will!

Alternative density tree plotting method in phytools (that also allows stochastic map trees)

A few of days ago, an R-sig-phylo subscriber asked if there was a way to represent the distribution of trees in a single plot in a manner “graphical style similar to the phytools::fancyTree(type = "phenogram95”)“.

In fact, such a function (densitTree) exists in phangorn; however since the subscriber (Giulio Dalla Riva at the University of Canterbury) added the challenge "bonus point if anybody proposes a solution that works with plotSimmap for painted trees”, I thought it might be worth demonstrating how this can be adapted to work with phytools plotting functions such as plotTree or plotSimmap.

I sent some code to the list, but it is fairly straightforward to write a custom function that does this visualization:

make.transparent<-function(color,alpha){
    RGB<-col2rgb(color)[,1]/255
    rgb(RGB[1],RGB[2],RGB[3],alpha)
}

density.tree<-function(trees,colors="blue",alpha=NULL,method="plotTree",fix.depth=FALSE,...){
    N<-length(trees)
    if(hasArg(use.edge.length)) use.edge.length<-list(...)$use.edge.length
    else use.edge.length<-TRUE
    if(!use.edge.length) trees<-lapply(trees,compute.brlen)
    if(!fix.depth){
        h<-sapply(trees,function(x) max(nodeHeights(x)))
        ii<-order(h,decreasing=TRUE)
        trees<-trees[ii]
        h<-h[ii]
    }
    if(is.null(alpha)) alpha<-1/N
    colors<-setNames(sapply(colors,make.transparent,alpha),names(colors))
    if(method=="plotTree") foo<-plotTree else foo<-plot
    foo(trees[[1]],color=colors,...)
    xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
    xlim[1]<-xlim[1]+0.03703704*diff(xlim)
    xlim[2]<-xlim[2]-0.03703704*diff(xlim)
    par(fg="transparent")
    for(i in 2:length(trees))
        foo(trees[[i]],
            tips=setNames(1:Ntip(trees[[1]]),trees[[1]]$tip.label),
            color=colors,add=TRUE,
            xlim=if(fix.depth) NULL else xlim-(h[1]-h[i]),...)
    par(fg="black")
}

First, let's use the function to do a basic visualization of topological and branch length uncertainty:

library(phytools)
trees<-read.nexus("Wayqecha.tree")
density.tree(trees,lwd=3,ftype="i",fsize=0.9)

We have all the options of plotTree & plotSimmap at our disposal here, and I find this type of visualization more readable with a non-standard node position which in phytools is called nodes="inner":

density.tree(trees,lwd=3,ftype="i",fsize=0.9,nodes="inner")

In this visualization, trees have variable depth (because the input trees were variable in total depth, but we can also turn this off and plot the trees with a standardized depth:

density.tree(trees,lwd=3,colors="black",ftype="i",fsize=0.9,nodes="inner",fix.depth=TRUE)

Finally, we can also turn off edge lengths entirely for this basic visualization by setting use.edge.length=FALSE:

density.tree(trees,lwd=3,ftype="i",fsize=0.9,nodes="inner",use.edge.length=FALSE)

As noted above, it is fairly straightforward to switch the plotting method used internally by this function to plotSimmap and thus plot stochastic map-style painted trees. For this example I will just use some simulated data:

Q<-matrix(c(-0.02,0.02,0.02,-0.02),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
x<-sim.history(trees[[1]],Q)$states

## Done simulation(s).

mtrees<-make.simmap(trees,x,message=FALSE)
mtrees

## 100 phylogenetic trees with mapped discrete characters

density.tree(mtrees,colors=setNames(c("red","blue"),sort(unique(x))),
    nodes="inner",ftype="i",lwd=5,method="plotSimmap")
add.simmap.legend(colors=setNames(c("red","blue"),sort(unique(x))),
    prompt=FALSE,x=par()$usr[1]+8,y=par()$usr[4]-1)

## or
density.tree(mtrees,colors=setNames(c("red","blue"),sort(unique(x))),
    nodes="inner",ftype="i",lwd=5,method="plotSimmap",fix.depth=TRUE)
add.simmap.legend(colors=setNames(c("red","blue"),sort(unique(x))),
    prompt=FALSE,x=par()$usr[1]+1,y=par()$usr[4]-1)

That's it really.

One feature that would be good to include is a preliminary calculation of the tip order in a consensus tree. Presently, the tip order is the order of the deepest tree - which is fine when the trees vary little in topology, as in this case, but will work poorly if that tree happens to be quite topologically different from the rest of the trees in the sample.