## Sunday, October 29, 2017

### New fix to optimization in `ratebytree` for comparing diversification rates

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.

#### 1 comment:

1. The handle was too replcia watches uk , and it was rolex replica ; so the founder of Xiaoshun personally remodeled the bag for his mother and named it “I AM NOT”, then decided to create a practical and design brand of the same name.