Graphing the likelihood surface of a fitted Mk discrete character evolution model using phytools

Here's a quick tutorial on graphing the likelihood surface for an Mk model using the phytools function fitMk.

Essentially, our workflow is that we first fit the model using fitMk. It's not actually important that we find the true MLE in this step because what we end up doing is pulling the likelihood function from this model object!

Next, we use this likelihood function to compute the log-likelihood for the values of the transition rate or rates, q, over our range of interest.

Finally, we'll select an appropriate plotting method and graph our surface.

For this demo I'll use data for feeding mode (a binary trait) modified from Revell & Collar (2009).

library(phytools)
data(sunfish.tree)
data(sunfish.data)
## extract discrete character (feeding mode)
fmode<-setNames(sunfish.data$feeding.mode,
    rownames(sunfish.data))

First, I'll compute a likelihood surface for a simple ER (equal-rates) model. This model has only one parameter, so this likelihood surface will have q on the horizontal axis, and log(L) on the vertical.

I'm going to define the range of values for q based on my pre-existing knowledge about the ML rates in this model. You should adjust to your own data accordingly.

## create dummy Q matrix of 1s
Q<-matrix(c(-1,1,1,-1),2,2,dimnames=list(levels(fmode),
    levels(fmode)))
## get likelihood function
fit<-fitMk(sunfish.tree,fmode,pi="fitzjohn")
lik<-fit$lik
## create vector
q<-seq(0.01,20,length.out=100)
## compute log-likelihoods
logL<-sapply(q,function(q) lik(Q*q))
## plot the surface
plot(q,logL,las=1,cex.axis=0.8,type="l",ylab="log(likelihood)",
    bty="n",las=1,col="grey",lwd=2)
## add the MLE back in
lines(rep(fit$rates,2),c(par()$usr[3],logLik(fit)),lty="dotted")
points(fit$rates,logLik(fit),pch=21,bg=palette()[4],cex=1.2)
legend("topright","MLE solution",pch=21,pt.bg=palette()[4],
    pt.cex=1.2,bg="white",box.col="white")

Second, let's repeat this but for an ARD (all-rates-different) model.

In this case, we have a binary trait, so this just involves estimating two different parameters, and then visualizing their joint surface.

logL<-matrix(NA,length(q),length(q),dimnames=list(q,q))
for(i in 1:nrow(logL)) for(j in 1:ncol(logL)) 
    logL[i,j]<-lik(matrix(c(-q[i],q[j],q[i],-q[j]),2,2))
contour(q,q,logL,nlevels=100,bty="n",col="grey",lwd=2)
title(xlab=expression(q["1,2"]),ylab=expression(q["2,1"]))
xy<-fitMk(sunfish.tree,fmode,model="ARD",pi="fitzjohn")$rates[2:1]
points(xy[1],xy[2],cex=1.5,pch=21,bg=palette()[4])
legend("topright","MLE solution",pch=21,pt.bg=palette()[4],pt.cex=1.2,
    bg="white",box.col="white")

Cool.

New co-phylogenetic method for phytools

For a day or two, I've been playing around with a new co-phylogenetic plotting method.

Here's what it looks like right now:

cotangleplot<-function(tr1,tr2,type=c("cladogram","phylogram"),
    use.edge.length=TRUE,tangle=c("both","tree1","tree2"),...){
    tr1<-untangle(tr1,"read.tree")
    tr2<-untangle(tr2,"read.tree")
    type<-type[1]
    tangle<-tangle[1]
    if(!use.edge.length){
        tr1<-compute.brlen(tr1)
        tr2<-compute.brlen(tr2)
    }
    if(hasArg(layout)) layout<-list(...)$layout
    else layout<-c(0.45,0.1,0.45)
    if(hasArg(nodes)) nodes<-list(...)$nodes
    else nodes<-"centered"
    if(hasArg(lty)) lty<-list(...)$lty
    else lty<-if(type=="phylogram") c("dotted","solid") else
        if(type=="cladogram") "solid"
    if(type=="phylogram") if(length(lty)==1) lty<-rep(lty,2)
    if(hasArg(lwd)) lwd<-list(...)$lwd
    else lwd<-2
    if(hasArg(cex)) cex<-list(...)$cex
    else cex<-1
    if(hasArg(color)) color<-list(...)$color
    else color<-palette()[4]
    if(tangle=="both"){
        capture.output(tmp<-cophylo(tr1,tr2))
        tips.tr1<-setNames(1:Ntip(tr1),tmp$trees[[1]]$tip.label)
        tips.tr2<-setNames(1:Ntip(tr2),tmp$trees[[2]]$tip.label)
        tips<-sort(rowMeans(cbind(tips.tr1,tips.tr2[names(tips.tr1)])))
        tips[]<-1:length(tips)
    } else if(tangle=="tree1"){
        tips<-setNames(1:Ntip(tr2),tr2$tip.label)
    } else if(tangle=="tree2"){
        tips<-setNames(1:Ntip(tr1),tr1$tip.label)
    }
    layout(matrix(c(1,2,3),1,3),widths=layout)
    plotTree(tr1,color="transparent",ftype="off",tips=tips,type=type,
        nodes=nodes)
    h<-par()$usr[2]
    pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    for(i in 1:Ntip(tr1)) lines(c(pp$xx[i],h),rep(pp$yy[i],2),
        lty="dotted")
    if(type=="phylogram"){
        for(i in 1:nrow(tr1$edge))
            lines(pp$xx[tr1$edge[i,]],rep(pp$yy[tr1$edge[i,2]],2),
                lwd=lwd,col=color,lty=lty[2])
        for(i in Ntip(tr1)+1:tr1$Nnode){
            dd<-Children(tr1,i)
            lines(rep(pp$xx[i],length(dd)),pp$yy[dd],lwd=lwd,
                col=color,lty=lty[1])   
        }
    } else if(type=="cladogram"){
        par(ljoin=2)
        for(i in 1:Ntip(tr1)){
            AA<-c(i,Ancestors(tr1,i))
            ii<-sapply(AA[1:(length(AA)-1)],function(x,y) 
                which(y==x), y=tr1$edge[,2])
            lines(c(pp$xx[tr1$edge[ii,2]],pp$xx[Ntip(tr1)+1]),
                c(pp$yy[tr1$edge[ii,2]],pp$yy[Ntip(tr1)+1]),
                lwd=lwd+2,
                col=if(par()$bg=="transparent") "white" else par()$bg,
                lty="solid")
            lines(c(pp$xx[tr1$edge[ii,2]],pp$xx[Ntip(tr1)+1]),
                c(pp$yy[tr1$edge[ii,2]],pp$yy[Ntip(tr1)+1]),
                lwd=lwd,col=color,lty=lty)
        }
    }   
    plot(NA,xlim=c(-1,1),ylim=pp$y.lim,axes=FALSE,xlab="",ylab="")
    text(rep(0,Ntip(tr1)),tips,gsub("_"," ",names(tips)),cex=cex,font=3)
    plotTree(tr2,color="transparent",ftype="off",direction="leftwards",
        tips=tips,type=type)
    h<-par()$usr[1]
    pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    for(i in 1:Ntip(tr2)) lines(c(pp$xx[i],h),rep(pp$yy[i],2),
        lty="dotted")
    if(type=="phylogram"){
        for(i in 1:nrow(tr2$edge))
            lines(pp$xx[tr2$edge[i,]],rep(pp$yy[tr2$edge[i,2]],2),
                lwd=lwd,col=color,lty=lty[2])
        for(i in Ntip(tr2)+1:tr2$Nnode){
            dd<-Children(tr2,i)
            lines(rep(pp$xx[i],length(dd)),sort(pp$yy[dd]),
                lwd=lwd,col=color,lty=lty[1])
        }
    } else if(type=="cladogram"){
        for(i in 1:Ntip(tr2)){
            AA<-c(i,Ancestors(tr2,i))
            ii<-sapply(AA[1:(length(AA)-1)],function(x,y) 
                which(y==x), y=tr2$edge[,2])
            lines(c(pp$xx[tr2$edge[ii,2]],pp$xx[Ntip(tr2)+1]),
                c(pp$yy[tr2$edge[ii,2]],pp$yy[Ntip(tr2)+1]),
                lwd=lwd+2,
                col=if(par()$bg=="transparent") "white" else par()$bg,
                lty="solid")
            lines(c(pp$xx[tr2$edge[ii,2]],pp$xx[Ntip(tr2)+1]),
                c(pp$yy[tr2$edge[ii,2]],pp$yy[Ntip(tr2)+1]),
                lwd=lwd,col=color,lty=lty)
        }
    }   
}

Let's give it a try using a phylogeny for salamanders (from Highton & Larson 1979) – and then a “perturbed” two different ways using random NNIs (using phangorn::rNNI. In practice, though, we might use this method to compare (for example) trees estimated using different methodologies or data for the same set of taxa.

set.seed(32)
library(phytools)
library(phangorn)
data(salamanders)
tree1<-force.ultrametric(rNNI(salamanders,3),
    message=FALSE)
tree2<-force.ultrametric(rNNI(salamanders,3),
    message=FALSE)

First let's see what the two trees look like plotted side-by-side.

par(mfrow=c(1,2))
plotTree(tree1)
plotTree(tree2)

Now let's compare them using our cotangleplot method.

cotangleplot(tree1,tree2,lwd=4,tangle="tree2")

This function gives us some different options. For instance, we can choose to tangle tree 1, tree 2, or both. Let's try a different tree of mammals from Garland et al. (1992), but set tangle="both".

data(mammal.tree)
par(bg="black",fg="white")
tree1<-force.ultrametric(rNNI(mammal.tree,3),
    message=FALSE)
tree2<-force.ultrametric(rNNI(mammal.tree,3),
    message=FALSE)
cotangleplot(tree1,tree2,lwd=4,tangle="tree2",use.edge.length=FALSE,
    layout=c(0.43,0.14,0.43),cex=0.9)

Finally, here's a square phylogram style, in which I make the plotted edges semi-transparent.

cotangleplot(tree1,tree2,lwd=3,tangle="tree1",type="phylogram",
    color=make.transparent(palette()[4],0.75),layout=c(0.43,0.14,0.43))

OK, it's not amazing, but it's not bad.