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

Make sure you only get the best replica parajumpers jackets to buy replica parajumpers kurtki damskie，fake parajumpers kurtki damskie puffer，discount parajumpers meska kurtka ，replica parajumpers jackets very similar to buying any other type...

ReplyDelete