Hi all,
This (ML lambda > mu) is a common belief that I think is mostly caused by the
fact that most implementations of Nee's equations constrain this.
library(diversitree)
set.seed(1)
trees <- replicate(40, rcoal(100), simplify=FALSE)
get.ml.bd <- function(phy, p=c(.1, 0))
coef(find.mle(make.bd(phy), p, method="subplex"))
ml.bd <- sapply(trees, get.ml.bd)
sum(ml.bd["lambda",] > ml.bd["mu",]) # 5/40
If you do not condition on survival of the two base lineages and the root
speciation event (which both Laser and diversitree do by default) then I
believe you'll always get ML lambda > mu (you will with the set of trees above,
anyway). A numerical demonstration that laser and diversitree's likelihoods
agree is below.
In an MCMC, when allowing lambda < mu, it is not mixing that is improved; you
are targeting a different distribution. When lambda > mu over most of the
posterior distribution, this will not make much difference, though.
However, there is a subtlety in allowing negative net diversification. If you
are working in {r, a} space and are using one-dimenensional parameter updates
(in an ML or MCMC search), you will never change from negative to positive net
diversification. This is because if r moves from a positive to a negative
value while a stays less than 1, this implies negative speciation and
extinction rates. If you work in {lambda, mu} space, this is not a problem.
All of this is really only an issue with coalescent-like trees, with a lot of
branching close to the present.
A small comment on the justification for the abs(...) that Emmanuel noted,
which is what allow calculation of the likelihoods with lambda < mu, is here:
http://www.zoology.ubc.ca/~fitzjohn/negative_r.pdf
(I wrote this two years ago and have not carefully checked it, so my apologies
if it is not totally clear.)
Cheers,
Rich
## For each tree, search for the ML point with Laser (so lambda>mu),
## then check that the diversitree's likelihood at this point agrees.
library(laser)
check <- function(phy) {
fit.laser <- bd(branching.times(phy))
## convert {r, a} -> {lambda, mu}
p <- with(fit.laser, c(r1/(1-a), a*r1/(1-a)))
lik <- make.bd(phy)
## Difference between diversitree and Laser at Laser's ML point:
lik(p) - fit.laser$LH
}
diff <- sapply(trees, check)
max(abs(diff)) # greatest difference is 8e-8, from roundoff error
On 2012-04-19, at 11:15 PM, Matt Pennell wrote:
> Hi all,
>
> Just a note on this topic, which Rich may be able to help clarify if he
> sees this. It can be shown that Rich's equation which uses log(abs(...)) is
> exactly equivalent to the regular version of the Nee equation. However, for
> an ultrametric tree, the ML of lambda will always be greater than mu. As I
> understand it, the advantage of the approach Diversitree uses is primarily
> in the context of using a mcmc chain to estimate the parameters
>
> myLikelihoodFunction <- make.bd(myTree, sampling.f = 1.0)
>
> p <- starting.point.bd(myTree)
>
> pars_init <- c(p[1], p[2])
>
> myprior <- make.prior.exponential(1 / (2 * (p[1] - p[2])))
>
> runmcmc <- mcmc(myLikelihoodFunction, x.init=pars_init, w=rep(1,2),
> prior=myprior)
>
> Allowing the chain to sample from parameter space where lambda < mu allows
> for better mixing of the mcmc.
>
> cheers,
> matt
>
> On Thu, Apr 19, 2012 at 9:29 PM, Emmanuel Paradis
> <emmanuel.para...@ird.fr>wrote:
>
>> Hi Simon,
>>
>> This is fixed and you can find the new version on ape's SVN. Another thing
>> I fixed: it could happen that the calculation of the Hessian-based SEs
>> fails: the error is now caught and the MLEs are returned with SEs set to NA
>> (the profile-likelihood CIs are still computed of course).
>>
>> Matt is right that the distribution of branching times is "weird" with a
>> coalescent tree, but an answer is still expected from birthdeath().
>>
>> Regarding the transformation on mu and lambda, the present version imposes
>> r (= lambda - mu) < 0 and/or a (= mu/lambda) > 1 because the likelihood
>> computation uses log(r) and log(1 - a). This also avoids most warnings now.
>> diversitree uses log(abs(... which, at first sight, is strange to me. But
>> there may be some justification for this. Rich may comment on it?
>>
>> Best,
>>
>> Emmanuel
>>
>> Simon Greenhill wrote on 20/04/2012 06:52:
>>
>> Morning all,
>>>
>>> Whenever I try to use ape's (v 3.0-2) birthdeath function e.g.
>>>
>>> birthdeath(rcoal(sample(10:**100, 1)))
>>>>
>>>
>>> I get the following error:
>>>
>>> Error in while (foo(up)< 0) up<- up + inc :
>>> missing value where TRUE/FALSE needed
>>>
>>> The warnings all appear to be this sort of thing:
>>>
>>> 1: In log(1 - a) : NaNs produced
>>> 2: In log(exp(r * x[2:N]) - a) : NaNs produced
>>> 3: In nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE) :
>>> NA/Inf replaced by maximum positive value
>>>
>>> and seem to be generated from this line:
>>> out<- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
>>>
>>>
>>> Is there a fix/workaround for this?
>>>
>>> Kind regards,
>>>
>>> Simon
>>>
>>> ______________________________**_________________
>>> R-sig-phylo mailing list
>>> R-sig-phylo@r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
>>>
>>>
>> --
>> Emmanuel Paradis
>> IRD, Jakarta, Indonesia
>> http://ape.mpl.ird.fr/
>>
>>
>> ______________________________**_________________
>> R-sig-phylo mailing list
>> R-sig-phylo@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
