Description Usage Arguments Details Author(s) Examples
View source: R/modelbissetd.R
Create a likelihood function for a BiSSE model where different chunks of time have different parameters. This code is experimental!
1 2 3 4 5 
tree 
An ultrametric bifurcating phylogenetic tree, in

states 
A vector of character states, each of which must be 0 or
1, or 
n.epoch 
Number of epochs. 1 corresponds to plain BiSSE, so this will generally be an integer at least 2. 
functions 
A named character vector of functions of time. See details. 
unresolved 
Unresolved clade information: see

sampling.f 
Vector of length 2 with the estimated proportion of
extant species in state 0 and 1 that are included in the phylogeny.
See 
nt.extra 
The number of species modelled in unresolved clades (this is in addition to the largest observed clade). 
strict 
The 
control 
List of control parameters for the ODE solver. See
details in 
truncate 
Logical, indicating if functions should be truncated to zero when negative (rather than causing an error). May be scalar (applying to all functions) or a vector (of length 6). 
spline.data 
List of data for splinebased time functions. See details 
.
This builds a BiSSE likelihood function where different regions of time (epochs) have different parameter sets. By default, all parameters are free to vary between epochs, so some constraining will probably be required to get reasonable answers.
For n
epochs, there are n1
time points; the first
n1
elements of the likelihood's parameter vector are these
points. These are measured from the present at time zero, with time
increasing towards the base of the tree. The rest of the parameter
vector are BiSSE parameters; the elements n:(n+6)
are for the
first epoch (closest to the present), elements (n+7):(n+13)
are
for the second epoch, and so on.
For make.bisse.t
, the funtions
is a vector of names of
functions of time.
For example, to have speciation rates be linear functions of
time, while the extinction and character change rates be constant with
respect to time, one can do
1 
The functions here must have t
as their first argument,
interpreted as time back from the present. Other possible functions
are "sigmoid.t", "stepf.t", "spline.t", "exp.t", and "spline.linear.t".
Unfortunately, documentation is still pending.
Richard G. FitzJohn
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  ## 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")
}
set.seed(4)
pars < c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
phy < tree.bisse(pars, max.t=30, x0=0)
## Suppose we want to see if diversification is different in the most
## recent 3 time units, compared with the rest of the tree (yes, this is
## a totally contrived example!):
plot(phy)
axisPhylo()
abline(v=max(branching.times(phy))  3, col="red", lty=3)
## For comparison, make a plain BiSSE likelihood function
lik.b < make.bisse(phy, phy$tip.state)
## Create the timedependent likelihood function. The final argument
## here is the number of 'epochs' that are allowed. Two epochs is one
## switch point.
lik.t < make.bisse.td(phy, phy$tip.state, 2)
## The switch point is the first argument. The remaining 12 parameters
## are the BiSSE parameters, with the first 6 being the most recent
## epoch.
argnames(lik.t)
pars.t < c(3, pars, pars)
names(pars.t) < argnames(lik.t)
## Calculations are identical to a reasonable tolerance:
lik.b(pars)  lik.t(pars.t)
## It will often be useful to constrain the time as a fixed quantity.
lik.t2 < constrain(lik.t, t.1 ~ 3)
## Parameter estimation under maximum likelihood. This is marked "don't
## run" because the timedependent fit takes a few minutes.
## Not run:
## Fit the BiSSE ML model
fit.b < find.mle(lik.b, pars)
## And fit the BiSSE/td model
fit.t < find.mle(lik.t2, pars.t[argnames(lik.t2)],
control=list(maxit=20000))
## Compare these two fits with a likelihood ratio test (lik.t2 is nested
## within lik.b)
anova(fit.b, td=fit.t)
## End(Not run)
## The time varying model (bisse.t) is more general, but substantially
## slower. Here, I will show that the two functions are equivalent for
## step function models. We'll constrain all the nonlambda parameters
## to be the same over a timeswitch at t=5. This leaves 8 parameters.
lik.td < make.bisse.td(phy, phy$tip.state, 2)
lik.td2 < constrain(lik.td, t.1 ~ 5,
mu0.2 ~ mu0.1, mu1.2 ~ mu1.1,
q01.2 ~ q01.1, q10.2 ~ q10.1)
lik.t < make.bisse.t(phy, phy$tip.state,
rep(c("stepf.t", "constant.t"), c(2, 4)))
lik.t2 < constrain(lik.t, lambda0.tc ~ 5, lambda1.tc ~ 5)
## Note that the argument names for these functions are different from
## one another. This reflects different ways that the functions will
## tend to be used, but is potentially confusing here.
argnames(lik.td2)
argnames(lik.t2)
## First, evaluate the functions with no time effect and check that they
## are the same as the base BiSSE model
p.td < c(pars, pars[1:2])
p.t < pars[c(1, 1, 2, 2, 3:6)]
## All agree:
lik.b(pars) # 159.7128
lik.td2(p.td) # 159.7128
lik.t2(p.t) # 159.7128
## In fact, the timevarying BiSSE will tend to be identical to plain
## BiSSE where the functions to not change:
lik.b(pars)  lik.t2(p.t)
## Slight numerical differences are typical for the timechunk BiSSE,
## because it forces the integration to be carried out more carefully
## around the switch point.
lik.b(pars)  lik.td2(p.td)
## Next, evaluate the functions with a time effect (5 time units ago,
## speciation rates were twice the contemporary rate)
p.td2 < c(pars, pars[1:2]*2)
p.t2 < c(pars[1], pars[1]*2, pars[2], pars[2]*2, pars[3:6])
## Huge drop in the likelihood (from 159.7128 to 172.7874)
lik.b(pars)
lik.td2(p.td2)
lik.t2(p.t2)
## The small difference remains between the two approaches, but they are
## basically the same.
lik.td2(p.td2)  lik.t2(p.t2)
## There is a small time cost to both timedependent methods,
## heavily paid for the timechunk case:
system.time(lik.b(pars))
system.time(lik.td2(p.td)) # 1.9x slower than plain BiSSE
system.time(lik.td2(p.td2)) # 1.9x slower than plain BiSSE
system.time(lik.t2(p.t)) # about the same speed
system.time(lik.t2(p.t2)) # about the same speed

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.