make.bisse: Binary State Speciation and Extinction Model

Description Usage Arguments Details Unresolved clade information ODE solver control Author(s) References See Also Examples

View source: R/model-bisse.R

Description

Prepare to run BiSSE (Binary State Speciation and Extinction) on a phylogenetic tree and character distribution. This function creates a likelihood function that can be used in maximum likelihood or Bayesian inference.

Usage

1
2
3
make.bisse(tree, states, unresolved=NULL, sampling.f=NULL, nt.extra=10,
           strict=TRUE, control=list())
starting.point.bisse(tree, q.div=5, yule=FALSE)

Arguments

tree

An ultrametric bifurcating phylogenetic tree, in ape “phylo” format.

states

A vector of character states, each of which must be 0 or 1, or NA if the state is unknown. This vector must have names that correspond to the tip labels in the phylogenetic tree (tree$tip.label). For tips corresponding to unresolved clades, the state should be NA.

unresolved

Unresolved clade information: see section below for structure.

sampling.f

Vector of length 2 with the estimated proportion of extant species in state 0 and 1 that are included in the phylogeny. A value of c(0.5, 0.75) means that half of species in state 0 and three quarters of species in state 1 are included in the phylogeny. By default all species are assumed to be known.

nt.extra

The number of species modelled in unresolved clades (this is in addition to the largest observed clade).

control

List of control parameters for the ODE solver. See details below.

strict

The states vector is always checked to make sure that the values are 0 and 1 only. If strict is TRUE (the default), then the additional check is made that every state is present. The likelihood models tend to be poorly behaved where states are missing.

q.div

Ratio of diversification rate to character change rate. Eventually this will be changed to allow for Mk2 to be used for estimating q parameters.

yule

Logical: should starting parameters be Yule estimates rather than birth-death estimates?

Details

make.bisse returns a function of class bisse. This function has argument list (and default values)

1
2
3
    f(pars, condition.surv=TRUE, root=ROOT.OBS, root.p=NULL,
      intermediates=FALSE)
  

The arguments are interpreted as

starting.point.bisse produces a heuristic starting point to start from, based on the character-independent birth-death model. You can probably do better than this; see the vignette, for example. bisse.starting.point is the same code, but deprecated in favour of starting.point.bisse - it will be removed in a future version.

Unresolved clade information

This must be a data.frame with at least the four columns

These columns may be in any order, and additional columns will be ignored. (Note that column names are case sensitive).

An alternative way of specifying unresolved clade information is to use the function make.clade.tree to construct a tree where tips that represent clades contain information about which species are contained within the clades. With a clade.tree, the unresolved object will be automatically constructed from the state information in states. (In this case, states must contain state information for the species contained within the unresolved clades.)

ODE solver control

The differential equations that define the BiSSE model are solved numerically using ODE solvers from the GSL library or deSolve's LSODA. The control argument to make.bisse controls the behaviour of the integrator. This is a list that may contain elements:

deSolve is the only supported backend on Windows.

Author(s)

Richard G. FitzJohn

References

FitzJohn R.G., Maddison W.P., and Otto S.P. 2009. Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies. Syst. Biol. 58:595-611.

Maddison W.P., Midford P.E., and Otto S.P. 2007. Estimating a binary character's effect on speciation and extinction. Syst. Biol. 56:701-710.

Nee S., May R.M., and Harvey P.H. 1994. The reconstructed evolutionary process. Philos. Trans. R. Soc. Lond. B Biol. Sci. 344:305-311.

See Also

constrain for making submodels, find.mle for ML parameter estimation, mcmc for MCMC integration, and make.bd for state-independent birth-death models.

The help pages for find.mle has further examples of ML searches on full and constrained BiSSE models.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
## Due to a change in sample() behaviour in newer R it is necessary to
## use an older algorithm to replicate the previous examples
if (getRversion() >= "3.6.0") {
  RNGkind(sample.kind = "Rounding")
}
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)

## Not run: 
## 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")

## End(Not run)

## 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)

diversitree documentation built on June 11, 2021, 5:17 p.m.