mbd_ml: Maximization of the loglikelihood under a multiple...
In Giappo/mbd: Multiple Birth Death Diversification

Description Usage Arguments Value Author(s) Examples

mbd_ml computes the maximum likelihood estimates of the parameters of a multiple birth-death diversification model for a given set of phylogenetic branching times. It also outputs the corresponding loglikelihood that can be used in model comparisons. Differently from mbd_ml it can account for three kind of events: sympatric (single) speciation, multiple (allopatric) speciation and extinction.

mbd_ml(
  loglik_function = mbd::mbd_loglik,
  brts,
  start_pars = c(0.5, 0.3, 0.5, 0.3),
  n_0 = 2,
  cond = 1,
  missnumspec = 0,
  lx = mbd::default_lx(brts = brts, missnumspec = missnumspec),
  ml_steps = 2,
  optim_ids = rep(TRUE, length(start_pars)),
  true_pars = start_pars,
  tips_interval = c(n_0 * (cond > 0), Inf),
  q_threshold = 0.001,
  verbose = TRUE,
  maxiter = 10000
)

`loglik_function`	the loglik function
`brts`	A set of branching times of a phylogeny.
`start_pars`	starting parameter values for the ml process.
`n_0`	the number of lineages at time equals zero.
`cond`	sets the conditioning cond = 0: no conditioning; cond = 1: conditioning on stem or crown age and non-extinction of the phylogeny;
`missnumspec`	The number of species that are in the clade, but missing in the phylogeny.
`lx`	it is the number of ODEs considered for the computation.
`ml_steps`	number of increasing precision steps to adopt in the ml routine. Each step will increment the absolute and relative tolerances of the optimizer as well as the number of equations used by the likelihood function.
`optim_ids`	ids of the parameters you want to optimize.
`true_pars`	true parameter values when running the ml process. The parameters that are not meant to be optimized specified as FALSE in `optim_ids` will be fixed by the homologous from this vector.
`tips_interval`	sets tips boundaries constrain on simulated dataset. It works only if cond == 1, otherwise it must be set to c(0, Inf).
`q_threshold`	adds a threshold for the evaluation of q. This is due because you never want `q` to actually be equal to zero or one.
`verbose`	choose if you want to print the output or not
`maxiter`	maximum number of subplex iterations

The output is a dataframe containing estimated parameters and maximum loglikelihood. The computed loglikelihood contains the factor q! m! / (q + m)! where q is the number of species in the phylogeny and m is the number of missing species, as explained in the supplementary material to Etienne et al. 2012.

Giovanni Laudanno

set.seed(2)
lambda <- 0.2 # sympatric speciation rate
mu <- 0.15 # extinction rate;
nu <- 2.0 # multiple allopatric speciation trigger rate
q <- 0.1 # single-lineage speciation probability
sim_pars <- c(lambda, mu, nu, q)
crown_age <- 1
cond <- 1
n_0 <- 2
sim <- mbd::mbd_sim(
 pars = sim_pars,
 n_0 = n_0, # Use a crown age
 age = crown_age,
 cond = cond # Condition on non-extinction
)
start_pars <- c(0.2, 0.15, 2, 0.15)
optim_ids <- c(FALSE, FALSE, FALSE, TRUE)
graphics::plot(sim$reconstructed_tree)
out <- mbd::mbd_ml(
  start_pars = start_pars,
  optim_ids = optim_ids,
  brts = sim$brts,
  true_pars = sim_pars,
  cond = cond,
  n_0 = n_0,
  verbose = FALSE
)