I just updated
(1,
2)
the *phytools* function `ratebytree`

for comparing the rates
of speciation & extinction between trees.

I noticed that (as I described
earlier
today for `fit.bd`

) occasionally during optimization of
the different diversification models, `ratebytree`

would return
non-finite values of the log-likelihood.

Since this almost always seems to occur far from the genuine ML values, the easiest thing to do was to merely add a check to verify that the log-likelihood was finite, and, if not, repeat optimization with different starting values for the parameters of the model until the optimization succeeds or until a certain number of iterations has been reached.

It seems to work. Here I'll demonstrate with the old function version -
then I'll load from source & show that it works. Obviously, the user should
just update *phytools* from GitHub rather than loading the function
from source!

```
library(phytools)
library(TreeSim)
```

First, a set of trees that fails under `model="equal-extinction"`

:

```
set.seed(304)
b<-(log(50)-log(2))/100
trees<-sim.bd.taxa.age(n=50,numbsim=3,lambda=b,mu=0,frac=1,age=100,T)
class(trees)<-"multiPhylo"
ratebytree(trees,model="equal-extinction")
```

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

```
## ML common diversification-rate model:
## b d k log(L)
## value 0.0385 0.0108 2 -200.9624
##
## ML multi diversification-rate model:
## b[1] b[2] b[3] d[1] d[2] d[3] k log(L)
## value 0.1527 0.3209 0.7647 0.6047 0.6047 0.6047 4 Inf
##
## Diversification model was "equal-extinction".
## Model fitting method was "nlminb".
##
## Likelihood ratio: Inf
## P-value (based on X^2): 0
```

```
## load from source
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/ratebytree.R")
ratebytree(trees,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
```

```
## ML common diversification-rate model:
## b d k log(L)
## value 0.0385 0.0108 2 -200.9624
##
## ML multi diversification-rate model:
## b[1] b[2] b[3] d[1] d[2] d[3] k log(L)
## value 0.0379 0.0383 0.0392 0.0107 0.0107 0.0107 4 -200.9466
##
## Diversification model was "equal-extinction".
## Model fitting method was "nlminb".
##
## Likelihood ratio: 0.0317
## P-value (based on X^2): 0.9843
```

Now a set that fails under `model="equal-speciation"`

:

```
rm(list=ls())
set.seed(9)
b<-(log(50)-log(2))/100
trees<-sim.bd.taxa.age(n=50,numbsim=3,lambda=b,mu=0,frac=1,age=100,T)
class(trees)<-"multiPhylo"
ratebytree(trees,model="equal-speciation")
```

```
## Warning in nlminb(c(init.b, rep(0, length(trees))), lik.eqlam, t = t, rho =
## rho, : NA/NaN function evaluation
```

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

```
## ML common diversification-rate model:
## b d k log(L)
## value 0.0332 0.0017 2 -204.3688
##
## ML multi diversification-rate model:
## b[1] b[2] b[3] d[1] d[2] d[3] k log(L)
## value 0.3515 0.3515 0.3515 0.9476 0 0 4 Inf
##
## Diversification model was "equal-speciation".
## Model fitting method was "nlminb".
##
## Likelihood ratio: Inf
## P-value (based on X^2): 0
```

```
## load from source
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/ratebytree.R")
ratebytree(trees,model="equal-speciation")
```

```
## Warning in nlminb(c(init.b, rep(0, length(trees))), lik.eqlam, t = t, rho =
## rho, : NA/NaN function evaluation
```

```
## Warning in nlminb(runif(n = 2, 0, 2) * c(init.b, rep(0, length(trees))), :
## NA/NaN function evaluation
```

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

```
## ML common diversification-rate model:
## b d k log(L)
## value 0.0332 0.0017 2 -204.3688
##
## ML multi diversification-rate model:
## b[1] b[2] b[3] d[1] d[2] d[3] k log(L)
## value 0.0334 0.0334 0.0334 0.0058 0 0 4 -204.1787
##
## Diversification model was "equal-speciation".
## Model fitting method was "nlminb".
##
## Likelihood ratio: 0.3803
## P-value (based on X^2): 0.8268
```

Since fixing `fit.bd`

I am not able to find failing cases for
`model="birth-death"`

, and I also could not find failing cases
for the common-rates model, but I added the fix anyway.

That's it.

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.