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)
```

## No comments:

## Post a Comment