Pybus & Harvey's γ through time

A phytools user today contacted me to ask about an evident discrepancy between an LTT plot (which on average seemed to be curved upward), and a calculated value of Pybus & Harvey's (2000) γ statistic (which was nonsignificant but negative). Visually, to me it seemed that this was most likely explained by a notable slightly negative curvature of the LTT plot towards the present day where there are many more internode distances. However, it also occurred to me that we can easily plot γ through time as well, and this might help us to see what parts of our plot are influencial in the measured value of γ.

Note that this assumes a pure-birth process of lineage accumulation through time. If we have extinction (which, in general, should result in positive γ), we need the 'pull of the present' to detect it, and by slicing of the end of our tree in computing γ we might find a value close to zero rather than positive!

Here are our functions:

gtt<-function(tree,n=100,...){
    if(hasArg(plot)) plot<-list(...)$plot
    else plot<-FALSE
    obj<-ltt(tree,plot=FALSE)
    t<-obj$times[which(obj$ltt==3)[1]]
    h<-max(nodeHeights(tree))
    x<-seq(t,h,by=(h-t)/(n-1))
    trees<-lapply(x,treeSlice,tree=tree,orientation="rootwards")
    gamma<-sapply(trees,function(x,plot){ 
        obj<-unlist(gammatest(ltt<-ltt(x,plot=FALSE)));
        if(plot) plot(ltt,xlim=c(0,h),ylim=c(1,Ntip(tree)),
            log.lineages=FALSE,log="y");
        Sys.sleep(0.01);
        obj},plot=plot)
    object<-list(t=x,gamma=gamma[1,],p=gamma[2,],tree=tree)
    class(object)<-"gtt"
    object
}
## plot method
plot.gtt<-function(x,...){
    args<-list(...)
    args$x<-x$t
    args$y<-x$gamma
    if(!is.null(args$show.tree)){ 
        show.tree<-args$show.tree
        args$show.tree<-NULL
    } else show.tree<-TRUE
    if(is.null(args$xlim)) args$xlim<-c(0,max(x$t))
    if(is.null(args$xlab)) args$xlab<-"time"
    if(is.null(args$ylab)) args$ylab<-expression(gamma)
    if(is.null(args$lwd)) args$lwd<-3
    if(is.null(args$type)) args$type<-"s"
    do.call(plot,args)
    if(show.tree) plotTree(x$tree,add=TRUE,ftype="off",mar=par()$mar,
        xlim=args$xlim,color=make.transparent("blue",0.1))
}
print.gtt<-function(x,...)
    cat("Object of class \"gtt\".\n\n")

Let's try it with a simulated tree:

library(phytools)
tree

## 
## Phylogenetic tree with 500 tips and 499 internal nodes.
## 
## Tip labels:
##  t106, t107, t97, t451, t452, t75, ...
## 
## Rooted; includes branch lengths.

object<-gtt(tree)
object

## Object of class "gtt".

plot(object)

Note that since our intervals are evenly spaced, rather than linked to events on the tree, the choice of a 'step' line plot is purely cosmetic and I could've also used straight lines:

plot(object,type="l")

The function gtt has an option plot=TRUE that will create an animation of the LTT plot for the tree. Something like this:

object<-gtt(tree,plot=TRUE)

The .gif was generated in R using ImageMagick as follows:

png(file="gtt-%03d.png",width=600,height=600)
object<-gtt(tree,plot=TRUE)
dev.off()
system("ImageMagick convert -delay 10 -loop 0 *.png gtt-anim.gif")
file.remove(list.files(pattern=".png"))

Type I error rates for variable-rate/process Mk discrete character evolution model

Earlier today I described a new method for testing hypotheses about heterogeneity in the rate (or process) of discrete character evolution across a tree.

Although mathematically quite different, this methods owes its intellectual provenance to an article by Brian O'Meara & colleagues published nearly a dozen years ago in 2006. The tactic here is basically the same, but applied to discrete rather than continuous traits.

I just pushed another small update that modifies the logLik methods of the "fitMk" and "fitmultiMk" object classes to facilitate comparison of models in which the rate (or process) is either homogeneous ("fitMk") or non-homogeneous ("fitmultiMk") across the edges of the tree.

In ensuring too that full credit goes where it is due, I should note that for computing the likelihood I adapted Felsenstein's pruning algorithm as implemented in Emmanuel Paradis's ape package.

Below, I'll examine the type I error of the method when used in hypothesis testing against a homogeneous rate/process model.

First, we need to simulate our 'regimes' on the trees. These could be specified arbitrarily, but in this case I will just simulate a binary character on a set of pure-birth phylogenies using sim.history:

library(phytools)
packageVersion("phytools")

## [1] '0.6.51'

trees<-pbtree(n=100,scale=1,nsim=200)
Q<-matrix(c(-0.5,0.5,0.5,-0.5),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
trees<-lapply(trees,sim.history,Q=Q,message=FALSE)
class(trees)<-c("multiSimmap","multiPhylo")

I generated 200 trees. Let's plot the first 100 of these:

par(mfrow=c(10,10))
nulo<-sapply(trees[1:100],plot,colors=setNames(c("blue","red"),
    letters[1:2]),ftype="off",lwd=1)

Next, we can simulate a discrete character on each of these trees, but in which the transition process is unrelated to the mapped regimes - that is, in which the null hypothesis of no difference in rate or process exists:

Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-0:1
X<-lapply(trees,function(x,Q) getStates(sim.history(x,Q,
    message=FALSE),"tips"),Q=Q)

Now let's fit our two models to each data vector & tree:

fits.single<-mapply(fitMk,tree=trees,x=X,MoreArgs=list(model="ER"),
    SIMPLIFY=FALSE)
fits.multi<-mapply(fitmultiMk,tree=trees,x=X,MoreArgs=list(model="ER"),
    SIMPLIFY=FALSE)

We can look at a single fitted model of each type:

fits.single[[1]]

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -0.912078  0.912078
## 1  0.912078 -0.912078
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -41.127165 
## 
## Optimization method used was "nlminb"

fits.multi[[1]]

## Object of class "fitmultiMk".
## 
## Fitted value of Q[a]:
##           0         1
## 0 -0.888472  0.888472
## 1  0.888472 -0.888472
## 
## Fitted value of Q[b]:
##           0         1
## 0 -1.032995  1.032995
## 1  1.032995 -1.032995
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -41.116092 
## 
## Optimization method used was "nlminb"

We can also conduct a likelihood-ratio test on each pair of fitted models. I will do that using lmtest::lrtest as follows:

library(lmtest)
suppressWarnings(LR.test<-mapply(lrtest,fits.single,fits.multi,
    SIMPLIFY=FALSE))
## for example
LR.test[[1]]

## Likelihood ratio test
## 
## Model 1: dots[[1L]][[1L]]
## Model 2: dots[[2L]][[1L]]
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   1 -41.127                     
## 2   2 -41.116  1 0.0221     0.8817

Let's pull out the P-values from these tests & plot them:

P<-sapply(LR.test,function(x) x[["Pr(>Chisq)"]][2])
obj<-hist(P,20,plot=FALSE)
plot(obj$mids,obj$counts/sum(obj$counts),type="h",
    lwd=18,col=make.transparent("blue",0.4),lend=1,
    xlab="P-value",ylab="relative frequency",ylim=c(0,0.2))
title(main="P-values for LR-test for data simulated under the null")
abline(h=0.05,col=make.transparent("red",0.5),lwd=1,
    lty="dotted")
text(0.96,0.05,"0.05",pos=3)

Neat. This is close to what we'd hope for.

Fitting a variable-process model of discrete character evolution on the tree using phytools

Now for something a little different.

Today, I have built a new method that fits a model of discrete character evolution in which the transition matrix Q varies among different parts of the tree.

These might be edges of clades specified arbitrarily by the user (for instance, using paintBranches or paintSubTree in phytools), or they could be regimes mapped onto the tree using the procedure of stochastic character mapping.

The way I did this was pretty simple. I just took the function that implements Felsenstein's famous pruning algorithm to compute the likelihood, but then I modified so that it could use a different Q for different edges. The only complication was that we might like our regime to change along an edge rather than merely at a node. To solve that, I used the phytools function map.to.singleton to convert our "simmap" object with singleton nodes and only a single regime per edge. Problem solved.

Note that this function, like fitMk, uses code for the pruning algorithm adapted from Emmanuel Paradis' ape package.

Let's try it.

First, our data:

library(phytools)
packageVersion("phytools")

## [1] '0.6.50'

plot(tree,ftype="off",colors=setNames(c("blue","red"),
    mapped.states(tree)),xlim=c(0,1.05*max(nodeHeights(tree))))
tiplabels(pie=to.matrix(x,c(0,1)),piecol=c("black","white"),
    cex=0.4,offset=0.01)

So the idea is simply that we will fit a model in which the rate of transition between 0 & 1 (show here as black & white) depends on the state (red or blue) mapped onto the edges of the tree. Right?

Let's fit our model:

fitERmulti<-fitmultiMk(tree,x,model="ER")
fitERmulti

## Object of class "fitmultiMk".
## 
## Fitted value of Q[a]:
##           0         1
## 0 -0.646581  0.646581
## 1  0.646581 -0.646581
## 
## Fitted value of Q[b]:
##           0         1
## 0 -10.14495  10.14495
## 1  10.14495 -10.14495
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -35.883658 
## 
## Optimization method used was "nlminb"

Or an "ARD" model that differs between parts of the tree:

fitARDmulti<-fitmultiMk(tree,x,model="ARD")
fitARDmulti

## Object of class "fitmultiMk".
## 
## Fitted value of Q[a]:
##           0         1
## 0 -2.218063  2.218063
## 1  0.604146 -0.604146
## 
## Fitted value of Q[b]:
##           0         1
## 0 -25.36092  25.36092
## 1  12.73165 -12.73165
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -34.515745 
## 
## Optimization method used was "nlminb"

Of course, we can compare this, if we'd like, to a model with but a single regime on the tree. For instance:

fitER<-fitMk(tree,x,model="ER")
fitER

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -1.651121  1.651121
## 1  1.651121 -1.651121
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -41.613061 
## 
## Optimization method used was "nlminb"

This suggests that our data justifies the greater model complexity of multiple regimes on the tree. That's good, because we simulated them that way!

tree<-pbtree(n=100,tip.label=LETTERS,scale=0.5)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
tree<-sim.history(tree,Q,anc="a")
sim.tree<-as.phylo(tree)
q<-setNames(c(1,10),letters[1:2])
sim.tree$edge.length<-colSums(t(tree$mapped.edge[,letters[1:2]])*q)
rownames(Q)<-colnames(Q)<-0:1
x<-as.factor(sim.history(sim.tree,Q)$states)

(As of yet there is not function to simulate multiple Mk models in different parts of the tree - so what I did above is to stretch the edge lengths of the tree by regime, and simulate under a constant regime on the stretched tree.)