R: Bootstrap likelihood ratio test of protracted birth-death...

pbd_LR {PBD}

R Documentation

Bootstrap likelihood ratio test of protracted birth-death model of diversification

Description

This function computes the maximum likelihood and the associated estimates of the parameters of a protracted birth-death model of diversification for a given set of phylogenetic branching times. It then performs a bootstrap likelihood ratio test of the protracted birth-death (PBD) model against the constant-rates (CR) birth-death model. Finally, it computes the power of this test.

Usage

pbd_LR(
  brts,
  initparsoptPBD,
  initparsoptCR,
  missnumspec,
  outputfilename = NULL,
  seed = 42,
  endmc = 1000,
  alpha = 0.05,
  plotit = TRUE,
  parsfunc = c(function(t,pars) {pars[1]},
  function(t,pars) {pars[2]},
  function(t,pars) {pars[3]},
  function(t,pars) {pars[4]}),
  cond = 1,
  btorph = 1,
  soc = 2,
  methode = 'lsoda',
  n_low = 0,
  n_up = 0,
  tol = c(1E-6,1E-6,1E-6),
  maxiter = 2000,
  optimmethod = 'subplex',
  verbose = FALSE
)

Arguments

`brts`	A set of branching times of a phylogeny, all positive
`initparsoptPBD`	The initial values of the parameters that must be optimized for the protracted birth-death (PBD) model: b, mu and lambda
`initparsoptCR`	The initial values of the parameters that must be optimized for the constant-rates (CR) model: b and mu
`missnumspec`	The number of species that are in the clade but missing in the phylogeny
`outputfilename`	The name (and location) of the file where the output will be saved. Default is no save.
`seed`	The seed for the pseudo random number generator for simulating the bootstrap data
`endmc`	The number of bootstraps
`alpha`	The significance level of the test
`plotit`	Boolean to plot results or not
`parsfunc`	Specifies functions how the rates depend on time, default functions are constant functions
`cond`	Conditioning: cond == 0 : conditioning on stem or crown age cond == 1 : conditioning on stem or crown age and non-extinction of the phylogeny cond == 2 : conditioning on stem or crown age and on the total number of extant taxa (including missing species) cond == 3 : conditioning on the total number of extant taxa (including missing species) Note: cond == 3 assumes a uniform prior on stem age, as is the standard in constant-rate birth-death models, see e.g. D. Aldous & L. Popovic 2004. Adv. Appl. Prob. 37: 1094-1115 and T. Stadler 2009. J. Theor. Biol. 261: 58-66.
`btorph`	Sets whether the likelihood is for the branching times (0) or the phylogeny (1)
`soc`	Sets whether stem or crown age should be used (1 or 2)
`methode`	The numerical method used to solve the master equation, such as 'lsoda' or 'ode45'.
`n_low`	Sets the lower bound of the number of species on which conditioning should be done when cond = 2. Set this to 0 when conditioning should be done on precisely the number of species (default)
`n_up`	Sets the upper bound of the number of species on which conditioning should be done when cond = 2. Set this to 0 when conditioning should be done on precisely the number of species (default)
`tol`	Sets the tolerances in the optimization. Consists of: reltolx = relative tolerance of parameter values in optimization reltolf = relative tolerance of function value in optimization abstolx = absolute tolerance of parameter values in optimization
`maxiter`	Sets the maximum number of iterations in the optimization
`optimmethod`	Method used in optimization of the likelihood. Current default is 'subplex'. Alternative is 'simplex'.
`verbose`	if TRUE, explanatory text will be shown

Details

The output is a list with 3 elements:

Value

`brtsCR`	a list of sets of branching times generated under the constant-rates model using the ML parameters under the CR model
`brtsDD`	a list of sets of branching times generated under the protracted birth-death model using the ML parameters under the PBD model
`out`	a dataframe with the parameter estimates and maximum likelihoods for protracted birth-death and constant-rates models `$model` - the model used to generate the data. 0 = unknown (for real data), 1 = CR, 2 = PBD `$mc` - the simulation number for each model `$b_CR` - speciation rate estimated under CR `$mu_CR` - extinction rate estimated under CR `$LL_CR` - maximum likelihood estimated under CR `$conv_CR` - convergence code for likelihood optimization; conv = 0 means convergence `$b_PBD1` - speciation-initation rate estimated under PBD for first set of initial values `$mu_PB1` - extinction rate estimated under DD for first set of initial values `$lambda_PB1` - speciation-completion rate estimated under PBD for first set of initial values `$LL_PBD1` - maximum likelihood estimated under DD for first set of initial values `$conv_PBD1` - convergence code for likelihood optimization for first set of initial values; conv = 0 means convergence `$b_PBD2` - speciation-initation rate estimated under PBD for second set of initial values `$mu_PB2` - extinction rate estimated under DD for second set of initial values `$lambda_PB2` - speciation-completion rate estimated under PBD for second set of initial values `$LL_PBD2` - maximum likelihood estimated under DD for second set of initial values `$conv_PBD2` - convergence code for likelihood optimization for second set of initial values; conv = 0 means convergence `$LR` - likelihood ratio between DD and CR
`pvalue`	p-value of the test
`LRalpha`	Likelihood ratio at the signifiance level alpha
`poweroftest`	power of the test for significance level alpha

Author(s)

Rampal S. Etienne

References

- Etienne, R.S. et al. 2016. Meth. Ecol. Evol. 7: 1092-1099, doi: 10.1111/2041-210X.12565