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.