## Simulate a tree and character distribution. This is on a birth-death
## tree, with high rates of character evolution and an asymmetry in the
## character transition rates.
pars <- c(.1, .1, .03, .03, .1, .2)
set.seed(3)
phy <- trees(pars, "bisse", max.taxa=25, max.t=Inf, x0=0)[[1]]
## Here is the 25 species tree with the true character history coded.
## Red is state '1', which has twice the character transition rate of
## black (state '0').
h <- history.from.sim.discrete(phy, 0:1)
plot(h, phy)
## Maximum likelihood parameter estimation:
p <- c(.1, .1) # initial parameter guess
## Not run:
# lik <- make.mk2(phy, phy$tip.state)
# fit.mk2 <- find.mle(lik, p)
# coef(fit.mk2) # q10 >> q01
# logLik(fit.mk2) # -10.9057
#
# ## This can also be done using the more general Mk-n.
# ## This uses an approximation for the likelihood calculations. make.mkn
# ## assumes that states are numbered 1, 2, ..., k, so 1 needs to be added
# ## to the states returned by trees.
# lik.mkn <- make.mkn(phy, phy$tip.state + 1, 2)
# fit.mkn <- find.mle(lik.mkn, p)
# fit.mkn[1:2]
#
# ## These are the same (except for the naming of arguments)
# all.equal(fit.mkn[-7], fit.mk2[-7], check.attr=FALSE, tolerance=1e-7)
#
# ## Equivalence to ape's \link{ace} function:
# model <- matrix(c(0, 2, 1, 0), 2)
# fit.ape <- ace(phy$tip.state, phy, "discrete", model=model, ip=p)
#
# ## To do the comparison, we need to rerun the diversitree version with
# ## the same root conditions as ape.
# fit.mk2 <- find.mle(lik, p, root=ROOT.GIVEN, root.p=c(1,1))
#
# ## These are the same to a reasonable degree of accuracy, too (the
# ## matrix exponentiation is slightly less accurate than the ODE
# ## solving approach. The make.mk2 version is exact)
# all.equal(fit.ape[c("rates", "loglik")], fit.mk2[1:2],
# check.attributes=FALSE, tolerance=1e-4)
#
# ## The ODE calculation method may be useful when there are a large
# ## number of possible states (say, over 20).
# lik.ode <- make.mkn(phy, phy$tip.state + 1, 2,
# control=list(method="ode"))
# fit.ode <- find.mle(lik.ode, p)
# fit.ode[1:2]
#
# all.equal(fit.ode[-7], fit.mkn[-7], tolerance=1e-7)
# ## End(Not run)
Run the code above in your browser using DataLab