make.musse.td {diversitree} | 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())
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 |
Please see make.bisse
for further details.
Richard G. FitzJohn
Here we will start with the tree and three-state 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 time-dependent 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)
## [1] "t.1" "lambda1.1" "lambda2.1" "lambda3.1" "mu1.1"
## [6] "mu2.1" "mu3.1" "q12.1" "q13.1" "q21.1"
## [11] "q23.1" "q31.1" "q32.1" "lambda1.2" "lambda2.2"
## [16] "lambda3.2" "mu1.2" "mu2.2" "mu3.2" "q12.2"
## [21] "q13.2" "q21.2" "q23.2" "q31.2" "q32.2"
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)
## [1] -3.822e-07
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 time-dependent fit takes a few minutes.
[This section is not run by default] 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)
[Ends not run section]