Comparing the rate of diversification among trees under a 'Yule' or pure-birth model

I just pushed another update to the function ratebytree for comparing diversification rates among trees - in this case to permit users to fit the 'Yule' or pure-birth model. In this case, two models are fit & compared: one in which each tree is permitted to have its own speciation rate & a second where all trees are constrained to the same rate.

Here is an example with two sets of trees - one (t1) simulated under the null hypothesis of pure-birth & no difference in rate among trees, and the second (t2) simulated with a 25% higher diversification rate in tree 3 compared to trees 1 & 2:

library(phytools)
## loading from source, although we can also update from GitHub
source("https://raw.githubusercontent.com/liamrevell/phytools/cd787d598150ac216fdbe87a60efb105e3138f0e/R/ratebytree.R")
print(t1,details=T)

## 3 phylogenetic trees
## tree 1 : 110 tips
## tree 2 : 79 tips
## tree 3 : 44 tips

par(mfrow=c(2,2))
plotTree(t1,ftype="off",lwd=1)

ratebytree(t1,model="Yule")

## ML common diversification-rate model:
##          b       d   k   log(L)
## value    0.037   0   1   -183.4468
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0378  0.0396  0.0315  0   0   0   3   -182.6673
## 
## Diversification model was "Yule".
## Model fitting method was "optimize".
## 
## Likelihood ratio: 1.5591 
## P-value (based on X^2): 0.4586

print(t2,details=T)

## 3 phylogenetic trees
## tree 1 : 73 tips
## tree 2 : 111 tips
## tree 3 : 235 tips

par(mfrow=c(2,2))
plotTree(t2,ftype="off",lwd=1)

ratebytree(t2,model="Yule")

## ML common diversification-rate model:
##          b       d   k   log(L)
## value    0.0427  0   1   -19.3546
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0355  0.0371  0.0494  0   0   0   3   -14.5664
## 
## Diversification model was "Yule".
## Model fitting method was "optimize".
## 
## Likelihood ratio: 9.5763 
## P-value (based on X^2): 0.0083

Neat. We can also compare these to alternative models, such as, for instance, our multi-rate birth-death model, using our AIC method. For instance:

yule<-ratebytree(t2,model="Yule")
bd<-ratebytree(t2,model="birth-death")

## Warning in nlminb(c(init.b, init.d), lik.bd, t = T, rho = rho, lower =
## rep(0, : NA/NaN function evaluation

## Warning in nlminb(c(init.b, init.d), lik.bd, t = T, rho = rho, lower =
## rep(0, : NA/NaN function evaluation

## Warning in nlminb(c(init.b, init.d), lik.bd, t = T, rho = rho, lower =
## rep(0, : NA/NaN function evaluation

AIC(yule,bd)

##                             AIC df
## common-rate            40.70918  1
## multi-rate:Yule        35.13287  3
## multi-rate:birth-death 40.72992  6

Although in this case our simplest model is just our Yule model and we have kind of ignored our birth-death null model (though we shouldn't have - I'll fix that later).

The trees were simulated using phytools as follows:

t1<-pbtree(b=0.039,t=100,nsim=3)
t2<-c(pbtree(b=0.039,t=100,nsim=2),
    as.multiPhylo(pbtree(b=0.049,t=100)))