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.
