pbd_LR() R function from [PBD]

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

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.


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:

Returns

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

References

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

Author(s)

Rampal S. Etienne

pbd_LR function