Two more models for comparing the diversification rates between reconstructed phylogenies

I just pushed (also here) a significant addition to the phytools function ratebytree for type="diversification" which allows the user to compare diversification rates (speciation & extinction) between phylogenies (e.g., 1, 2).

This update now permits us to fit two additional models: one (model="equal-extinction") in which speciation is permitted to vary among trees, but not extinction; and a second (model="equal-speciation") in which extinction may vary, but not speciation.

In the following, I'll demonstrate this with two sets of simulated trees as follows:

library(phytools)
packageVersion("phytools") ## from GitHub

## [1] '0.6.36'

print(tt1,details=TRUE)

## 3 phylogenetic trees
## tree 1 : 500 tips
## tree 2 : 500 tips
## tree 3 : 500 tips

print(tt2,details=TRUE)

## 3 phylogenetic trees
## tree 1 : 500 tips
## tree 2 : 500 tips
## tree 3 : 500 tips

First, let's visualize the branching times in our first set of trees, tt1:

obj<-ltt(tt1,plot=FALSE)
xlim<-c(max(sapply(obj,function(x) max(x$times))),0)
ylim<-c(1,max(sapply(obj,function(x) max(x$ltt))))
plot(max(obj[[1]]$times)-obj[[1]]$times,obj[[1]]$ltt,xlim=xlim,
    ylim=ylim,log="y",type="s",xlab="time before present",
    ylab="lineages",col="red",lwd=2,lend=1)
lines(max(obj[[2]]$times)-obj[[2]]$times,obj[[2]]$ltt,
    col="orange",lwd=2,lend=1,type="s")
lines(max(obj[[3]]$times)-obj[[3]]$times,obj[[3]]$ltt,
    col="blue",lwd=2,lend=1,type="s")
text(max(obj[[1]]$times),1.1,"tree 1",pos=4)
text(max(obj[[2]]$times),1.1,"tree 2",pos=4)
text(max(obj[[3]]$times),1.1,"tree 3",pos=4)
title(main="Equal speciation / variable extinction")

Now let's fit our different models. In total, we have four: a model in which speciation & extinction rates are equal between trees, a second in which speciation rates vary but extinction rates are equal between phylogenies, a third in which the converse is true, and, finally, a fourth in which birth & death rates are both allowed to vary among phylogenies.

eqex<-ratebytree(tt1,model="equal-extinction")

## Warning in nlminb(c(init.b, 0), lik.eqmu, t = t, rho = rho, trace =
## trace, : NA/NaN function evaluation

## Warning in nlminb(c(init.b, init.d), lik.onerate, t = t, rho = rho, model =
## model, : NA/NaN function evaluation

eqex

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0491  0.0104  2   1652.323
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0524  0.0486  0.0465  0.0101  0.0101  0.0101  4   1654.3604
## 
## Diversification model was "equal-extinction".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 4.0748 
## P-value (based on X^2): 0.1304

eqsp<-ratebytree(tt1,model="equal-speciation")

## Warning in nlminb(c(init.b, init.d), lik.onerate, t = t, rho = rho, model =
## model, : NA/NaN function evaluation

eqsp

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0491  0.0104  2   1652.323
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0486  0.0486  0.0486  7e-04   0.0097  0.0163  4   1657.5991
## 
## Diversification model was "equal-speciation".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 10.5522 
## P-value (based on X^2): 0.0051

bd<-ratebytree(tt1)

## 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

bd

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0491  0.0104  2   1652.323
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0471  0.0475  0.0546  0   0.0081  0.0243  6   1658.7793
## 
## Diversification model was "birth-death".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 12.9127 
## P-value (based on X^2): 0.0117

AIC(eqex,eqsp,bd)

##                                   AIC df
## common-rate                 -3300.646  2
## multi-rate:equal-extinction -3300.721  4
## multi-rate:equal-speciation -3307.198  4
## multi-rate:birth-death      -3305.559  6

aic.w(setNames(AIC(eqex,eqsp,bd)[,1],
    c("common-rate","equal-extinction","equal-speciation",
    "birth-death")))

##      common-rate equal-extinction equal-speciation      birth-death 
##       0.02489281       0.02584185       0.65896525       0.29030009

This shows us that (by a small margin) the favored model is one in which all three trees share a common speciation rate - but differ in their fitted extinction rates.

Now let's try again with our second set of phylogenies. First plotting them:

print(tt2,details=TRUE)

## 3 phylogenetic trees
## tree 1 : 500 tips
## tree 2 : 500 tips
## tree 3 : 500 tips

obj<-ltt(tt2,plot=FALSE)
xlim<-c(max(sapply(obj,function(x) max(x$times))),0)
ylim<-c(1,max(sapply(obj,function(x) max(x$ltt))))
plot(max(obj[[1]]$times)-obj[[1]]$times,obj[[1]]$ltt,xlim=xlim,
    ylim=ylim,log="y",type="s",xlab="time before present",
    ylab="lineages",col="red",lwd=2,lend=1)
lines(max(obj[[2]]$times)-obj[[2]]$times,obj[[2]]$ltt,
    col="orange",lwd=2,lend=1,type="s")
lines(max(obj[[3]]$times)-obj[[3]]$times,obj[[3]]$ltt,
    col="blue",lwd=2,lend=1,type="s")
text(max(obj[[1]]$times),1.1,"tree 1",pos=4)
text(max(obj[[2]]$times),1.1,"tree 2",pos=4)
text(max(obj[[3]]$times),1.1,"tree 3",pos=4)
title(main="Equal extinction / variable speciation")

Now once again we can fit our four different models & compare them:

eqex<-ratebytree(tt2,model="equal-extinction")

## Warning in nlminb(c(init.b, init.d), lik.onerate, t = t, rho = rho, model =
## model, : NA/NaN function evaluation

eqex

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0612  0.0148  2   1951.3546
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0514  0.0599  0.0701  0.0115  0.0115  0.0115  4   1964.3627
## 
## Diversification model was "equal-extinction".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 26.0162 
## P-value (based on X^2): 0

eqsp<-ratebytree(tt2,model="equal-speciation")

## Warning in nlminb(c(init.b, init.d), lik.onerate, t = t, rho = rho, model =
## model, : NA/NaN function evaluation

eqsp

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0612  0.0148  2   1951.3546
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0604  0.0604  0.0604  0.0229  0.012   0   4   1958.951
## 
## Diversification model was "equal-speciation".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 15.1928 
## P-value (based on X^2): 5e-04

bd<-ratebytree(tt2)

## Warning in nlminb(c(init.b, init.d), lik.onerate, t = t, rho = rho, model =
## model, : NA/NaN function evaluation

bd

## ML common diversification-rate model:
##          b       d       k   log(L)
## value    0.0612  0.0148  2   1951.3546
## 
## ML multi diversification-rate model:
##          b[1]    b[2]    b[3]    d[1]    d[2]    d[3]    k   log(L)
## value    0.0489  0.0595  0.0756  0.0068  0.0106  0.0216  6   1965.0754
## 
## Diversification model was "birth-death".
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 27.4415 
## P-value (based on X^2): 0

AIC(eqex,eqsp,bd)

##                                   AIC df
## common-rate                 -3898.709  2
## multi-rate:equal-extinction -3920.725  4
## multi-rate:equal-speciation -3909.902  4
## multi-rate:birth-death      -3918.151  6

aic.w(setNames(AIC(eqex,eqsp,bd)[,1],
    c("common-rate","equal-extinction","equal-speciation",
    "birth-death")))

##      common-rate equal-extinction equal-speciation      birth-death 
##       0.00001294       0.78095308       0.00348628       0.21554770

This shows highest support for the 'equal-extinction' variable speciation model.

Cool.

I simulated the trees above using phytools::pbtree conditioning on both time and the total number of taxa. This is done via rejection sampling & is very slow, so in reality I would probably recommend using the CRAN package TreeSim by Tanja Stadler which is evidently much faster.

library(phytools)
## simulate with constant speciation, variable extinction
tt1<-c(pbtree(b=0.05,t=110.4,n=500,method="direct"),
    pbtree(b=0.05,d=0.01,t=138.0,n=500,extant.only=T),
    pbtree(b=0.05,d=0.02,t=184.0,n=500,extant.only=T))
print(tt1,details=T)
## simulate with constant extinction, variable extinction
tt2<-c(pbtree(b=0.05,d=0.01,t=138.0,n=500,extant.only=T),
    pbtree(b=0.06,d=0.01,t=110.4,n=500,extant.only=T),
    pbtree(b=0.07,d=0.01,t=92.0,n=500,extant.only=T))
print(tt2,details=T)