View source: R/modelmussetd.R
make.musse.td  R Documentation 
Create a likelihood function for a MuSSE model where different chunks of time have different parameters. This code is experimental!
make.musse.td(tree, states, k, n.epoch, sampling.f=NULL,
strict=TRUE, control=list())
make.musse.t(tree, states, k, functions, sampling.f=NULL,
strict=TRUE, control=list(), truncate=FALSE, spline.data=NULL)
tree 
An ultrametric bifurcating phylogenetic tree, in

states 
A vector of character states, each of which must be an
integer between 1 and 
k 
The number of states. 
n.epoch 
Number of epochs. 1 corresponds to plain MuSSE, so this will generally be an integer at least 2. 
functions 
A named list of functions of time. See details. 
sampling.f 
Vector of length 
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 same length as the functions vector). 
spline.data 
List of data for splinebased time functions. See details 
.
Please see make.bisse.t
for further details.
Richard G. FitzJohn
## Here we will start with the tree and threestate character set from
## the make.musse example. This is a poorly contrived example.
pars < c(.1, .15, .2, # lambda 1, 2, 3
.03, .045, .06, # mu 1, 2, 3
.05, 0, # q12, q13
.05, .05, # q21, q23
0, .05) # q31, q32
set.seed(2)
phy < tree.musse(pars, 30, x0=1)
## 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 MuSSE likelihood function
lik.m < make.musse(phy, phy$tip.state, 3)
## 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.musse.td(phy, phy$tip.state, 3, 2)
## The switch point is the first argument. The remaining 12 parameters
## are the MuSSE 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.m(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 MuSSE ML model
fit.m < find.mle(lik.m, pars)
## And fit the MuSSE/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.m)
anova(fit.m, td=fit.t)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.