pbd_LR: 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

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

Details

The output is a list with 3 elements:

References