pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
set.seed(4)
phy <- tree.bisse(pars, max.t=30, x0=0)
## Here is the 52 species tree with the true character history coded.
## Red is state '1', which has twice the speciation rate of black (state
## '0').
h <- history.from.sim.discrete(phy, 0:1)
plot(h, phy)
lik <- make.bisse(phy, phy$tip.state)
lik(pars) # -159.71
## Heuristic guess at a starting point, based on the constant-rate
## birth-death model (uses \link{make.bd}).
p <- starting.point.bisse(phy)
## Start an ML search from this point. This takes some time (~7s)
fit <- find.mle(lik, p, method="subplex")
logLik(fit) # -158.6875
## The estimated parameters aren't too far away from the real ones, even
## with such a small tree
rbind(real=pars,
estimated=round(coef(fit), 2))
## Test a constrained model where the speciation rates are set equal
## (takes ~4s).
lik.l <- constrain(lik, lambda1 ~ lambda0)
fit.l <- find.mle(lik.l, p[-1], method="subplex")
logLik(fit.l) # -158.7357
## Despite the difference in the estimated parameters, there is no
## statistical support for this difference:
anova(fit, equal.lambda=fit.l)
## Run an MCMC. Because we are fitting six parameters to a tree with
## only 50 species, priors will be needed. I will use an exponential
## prior with rate 1/(2r), where r is the character independent
## diversificiation rate:
prior <- make.prior.exponential(1 / (2 * (p[1] - p[3])))
## This takes quite a while to run, so is not run by default
tmp <- mcmc(lik, fit$par, nsteps=100, prior=prior, w=.1, print.every=0)
w <- diff(sapply(tmp[2:7], range))
samples <- mcmc(lik, fit$par, nsteps=1000, prior=prior, w=w,
print.every=100)
## See \link{profiles.plot} for more information on plotting these
## profiles.
col <- c("blue", "red")
profiles.plot(samples[c("lambda0", "lambda1")], col.line=col, las=1,
xlab="Speciation rate", legend="topright")
## BiSSE reduces to the birth-death model and Mk2 when diversification
## is state independent (i.e., lambda0 ~ lambda1 and mu0 ~ mu1).
lik.mk2 <- make.mk2(phy, phy$tip.state)
lik.bd <- make.bd(phy)
## 1. BiSSE / Birth-Death
## Set the q01 and q10 parameters to arbitrary numbers (need not be
## symmetric), and constrain the lambdas and mus to be the same for each
## state. The likelihood function now has just two parameters and
## will be proprtional to Nee's birth-death based likelihood:
lik.bisse.bd <- constrain(lik,
lambda1 ~ lambda0, mu1 ~ mu0,
q01 ~ .01, q10 ~ .02)
pars <- c(.1, .03)
## These differ by -167.3861 for both parameter sets:
lik.bisse.bd(pars) - lik.bd(pars)
lik.bisse.bd(2*pars) - lik.bd(2*pars)
## 2. BiSSE / Mk2
## Same idea as above: set all diversification parameters to arbitrary
## values (but symmetric this time):
lik.bisse.mk2 <- constrain(lik,
lambda0 ~ .1, lambda1 ~ .1,
mu0 ~ .03, mu1 ~ .03)
## Differ by -150.4740 for both parameter sets.
lik.bisse.mk2(pars) - lik.mk2(pars)
lik.bisse.mk2(2*pars) - lik.mk2(2*pars)
## 3. Sampled BiSSE / Birth-Death
## Pretend that the tree is only .6 sampled:
lik.bd2 <- make.bd(phy, sampling.f=.6)
lik.bisse2 <- make.bisse(phy, phy$tip.state, sampling.f=c(.6, .6))
lik.bisse2.bd <- constrain(lik.bisse2,
lambda1 ~ lambda0, mu1 ~ mu0,
q01 ~ .01, q10 ~ .01)
## Difference of -167.2876
lik.bisse2.bd(pars) - lik.bd2(pars)
lik.bisse2.bd(2*pars) - lik.bd2(2*pars)
## 4. Unresolved clade BiSSE / Birth-Death
set.seed(1)
unresolved <- data.frame(tip.label=I(sample(phy$tip.label, 5)),
Nc=sample(2:10, 5), n0=0, n1=0)
unresolved.bd <- structure(unresolved$Nc, names=unresolved$tip.label)
lik.bisse3 <- make.bisse(phy, phy$tip.state, unresolved)
lik.bisse3.bd <- constrain(lik.bisse3,
lambda1 ~ lambda0, mu1 ~ mu0,
q01 ~ .01, q10 ~ .01)
lik.bd3 <- make.bd(phy, unresolved=unresolved.bd)
## Difference of -167.1523
lik.bisse3.bd(pars) - lik.bd3(pars)
lik.bisse3.bd(pars*2) - lik.bd3(pars*2)
## This only works when diversitree was built with CVODES support
## (currently not on Windows). This support is still experimental.
## Different backends:
lik.1 <- make.bisse(phy, phy$tip.state, control=list(backend="deSolve"))
lik.2 <- make.bisse(phy, phy$tip.state, control=list(backend="cvodes"))
lik.3 <- make.bisse(phy, phy$tip.state, control=list(backend="CVODES"))
## All the calculations are very similar
pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
(ll.1 <- lik.1(pars))
(ll.2 <- lik.2(pars))
(ll.3 <- lik.3(pars))
diff(range(c(ll.1, ll.2, ll.3)))
## However, the CVODES version is quite a bit faster than the other two
## (cvodes is often a bit slower than deSolve).
## Timings are on a 2008 Mac Pro
system.time(fit.1 <- find.mle(lik.1, pars)) # 4.5s
system.time(fit.2 <- find.mle(lik.2, pars)) # 6.0s
system.time(fit.3 <- find.mle(lik.3, pars)) # 1.6s
## All found the same ML point:
diff(range(fit.1$lnLik, fit.2$lnLik, fit.3$lnLik))Run the code above in your browser using DataLab