Skip to content

Commit

Permalink
Simplify examples
Browse files Browse the repository at this point in the history
  • Loading branch information
richfitz committed Dec 13, 2023
1 parent a7fae7b commit f891037
Showing 1 changed file with 0 additions and 99 deletions.
99 changes: 0 additions & 99 deletions man/make.bisse.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -197,105 +197,6 @@ 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)

\dontrun{
## 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
unresolved <- data.frame(tip.label=I(c("sp25", "sp30", "sp40", "sp56", "sp20")),
Nc =c(10, 9, 6, 5, 2),
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)
}

\references{
Expand Down

0 comments on commit f891037

Please sign in to comment.