Package 'PBD' reference manual

Title:	Protracted Birth-Death Model of Diversification
Description:	Conducts maximum likelihood analysis and simulation of the protracted birth-death model of diversification. See Etienne, R.S. & J. Rosindell 2012 <doi:10.1093/sysbio/syr091>; Lambert, A., H. Morlon & R.S. Etienne 2014, <doi:10.1007/s00285-014-0767-x>; Etienne, R.S., H. Morlon & A. Lambert 2014, <doi:10.1111/evo.12433>.
Authors:	Rampal S. Etienne
Maintainer:	Rampal S. Etienne <[email protected]>
License:	GPL-2
Version:	1.5
Built:	2025-02-11 05:26:12 UTC
Source:	https://github.com/rsetienne/pbd

Protracted Birth-Death Model of Diversification

Description

This package computes the (maximum) likelihood of the protracted speciation model for a given set of branching times This package is a likelihood-based statistical package to estimate parameters under the protracted speciation model.

First version: 0.8
New in version 0.9
- Bug fix for stem age
New in version 0.91
- Reports loglik = -Inf on an error in the deSolve package (function ode)
New in version 0.92
- Correcting order of parameters of pbd_sim
New in version 0.93
- pbd_sim produces a tree, a matrix containing all events in the simulation, and a tree with one sample per species.
New in version 1.0
- Conditioning is also possible on a range of values of the number of species.
New in version 1.1
- Simulation of the protracted speciation tree has more features.
New in version 1.2
- Optimization can make use of subplex (default) and simplex (older versions).
New in version 1.3
- Contains a function to carry out a bootstrap likelihood ratio test.
- Vignette and test added.
- Reports an error if exteq = TRUE and initparsopt contains 4 parameters.
- Option to limit a simulation to a certain maximum number of species; if exceeded, the simulation is ignored.
New in version 1.4:
- Includes all special cases in pbd_durspec_mean
- Fixes a bug in conditioning on a range of values of the number of species
New in version 1.5:
- Fixes bug with missing species in dd_LR

Details

pbd_loglik computes the likelihood of the protracted birth-death model of diversification, given a set of parameters and a data set of phylogenetic branching times.

pbd_ML finds the parameters that maximizes the likelihood computed by pbd_loglik.

pbd_bootstrap performs a maximum likelihood analysis and simulates with the maximum likelihood parameters. The ML parameters of the simulated data sets are then estimated, providing an uncertainty distribution for the original ML estimate on the original data.

Author(s)

Rampal S. Etienne Maintainer: Rampal S. Etienne <[email protected]>

References

- Etienne, R.S. & J. Rosindell 2012. Systematic Biology 61: 204-213.
- Lambert, A., H. Morlon & R.S. Etienne 2014. Journal of Mathmematical Biology 70: 367-397. doi:10.1007/s00285-014-0767-x
- Etienne, R.S., H. Morlon & A. Lambert 2014. Evolution 68: 2430-2440, doi:10.1111/evo.12433
. - Etienne, R.S., A.L. Pigot & A.B. Phillimore 2016. Methods in Ecology & Evolution 7: 1092-1099, doi: 10.1111/2041-210X.12565

Bootstrap analysis under protracted birth-death model of diversification

Description

Likelihood maximization for protracted birth-death model of diversification followed by simulations of the model using the maximum likelihood parameter estimates to compute an estimate of the error in these estimates and to assess the goodness-of-fit of the model by comparing maximum likelihoods of the simulated data sets to the maximum likelihood of the real data set.

Usage

pbd_bootstrap(
  brts,
  initparsopt = c(0.2, 0.1, 1),
  idparsopt = 1:length(initparsopt),
  idparsfix = NULL,
  parsfix = NULL,
  exteq = (length(initparsopt) < 4),
  parsfunc = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  missnumspec = 0,
  cond = 1,
  btorph = 0,
  soc = 2,
  plotltt = 1,
  methode = "lsoda",
  n_low = 0,
  n_up = 0,
  tol = c(1e-04, 1e-04, 1e-06),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  endmc = 100,
  seed = 42,
  optimmethod = "subplex",
  num_cycles = 1,
  verbose = FALSE
)
pbd_bootstrap(
  brts,
  initparsopt = c(0.2, 0.1, 1),
  idparsopt = 1:length(initparsopt),
  idparsfix = NULL,
  parsfix = NULL,
  exteq = (length(initparsopt) < 4),
  parsfunc = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  missnumspec = 0,
  cond = 1,
  btorph = 0,
  soc = 2,
  plotltt = 1,
  methode = "lsoda",
  n_low = 0,
  n_up = 0,
  tol = c(1e-04, 1e-04, 1e-06),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  endmc = 100,
  seed = 42,
  optimmethod = "subplex",
  num_cycles = 1,
  verbose = FALSE
)

Arguments

`brts`	A set of branching times of a phylogeny, all positive
`initparsopt`	The initial values of the parameters that must be optimized
`idparsopt`	The ids of the parameters that must be optimized, e.g. 1:4 for all parameters. The ids are defined as follows: id == 1 corresponds to b (speciation-initiation rate) id == 2 corresponds to mu_1 (extinction rate of good species) id == 3 corresponds to la_1 (speciation-completion rate) id == 4 corresponds to mu_2 (extinction rate of incipient species)
`idparsfix`	The ids of the parameters that should not be optimized, e.g. c(2,4) if mu_1 and mu_2 should not be optimized, but only b and la_1. In that case idparsopt must be c(1,3).
`parsfix`	The values of the parameters that should not be optimized
`exteq`	Sets whether incipient species have the same (1) or different (0) extinction rate as good species. If exteq = 0, then idparsfix and idparsopt should together have all parameters 1:4
`parsfunc`	Specifies functions how the rates depend on time, default functions are constant functions
`missnumspec`	The number of species that are in the clade but missing in the phylogeny
`cond`	Conditioning: cond == 0 : conditioning on stem or crown age cond == 1 : conditioning on stem or crown age and non-extinction of the phylogeny
`btorph`	Sets whether the likelihood is for the branching times (0) or the phylogeny (1)
`soc`	Sets whether the first element of the branching times is the stem (1) or the crown (2) age
`plotltt`	Sets whether the lineage-through-time plot should be plotted (1) or not (0)
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.
`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
`endmc`	Sets the number of simulations for the bootstrap
`seed`	Sets the seed for the simulations of the bootstrap
`optimmethod`	Method used in optimization of the likelihood. Current default is 'subplex'. Alternative is 'simplex' (default of previous versions)
`num_cycles`	Number of cycles of the optimization (default is 1).
`verbose`	if TRUE, explanatory text will be shown

Value

A list of three dataframes. The first dataframe contains the maximum likelihood results of the real data set, the second contains the simulated trees, and the third dataframe, with number of rows equal to endmc, contain the maximum likelihood results for the simulated data. The columns of both frames contains the following elements for each simulated data set:

`ntips`	gives the number of tips
`b`	gives the maximum likelihood estimate of b
`mu_1`	gives the maximum likelihood estimate of mu_1
`la_1`	gives the maximum likelihood estimate of la_1
`mu_2`	gives the maximum likelihood estimate of mu_2
`loglik`	gives the maximum loglikelihood
`df`	gives the number of estimated parameters, i.e. degrees of feedom
`conv`	gives a message on convergence of optimization; conv = 0 means convergence
`exp_durspec`	gives the expected duration of speciation
`median_durspec`	gives the median duration of speciation

Author(s)

Rampal S. Etienne

Node depth probbaility density for protracted birth-death model of diversification

Description

pbd_brts_density computes the probability density of node depths under the protracted speciation model given a set of parameters

Usage

pbd_brts_density(
  pars1,
  pars1f = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  methode = "lsoda",
  brts
)
pbd_brts_density(
  pars1,
  pars1f = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  methode = "lsoda",
  brts
)

Arguments

`pars1`	Vector of parameters: `pars1[1]` corresponds to b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `pars1[2]` corresponds to mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `pars1[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars1[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species When rates depend on time this time dependence should be specified in pars1f and pars1 then becomes the parameters used in pars1f
`pars1f`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `pars1f[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `pars1f[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `pars1f[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars1f[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`methode`	sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.
`brts`	A set of branching times of a phylogeny, all positive, for which the density must be computed

Value

The probability density for all branching times

Author(s)

Rampal S. Etienne

Examples

 pbd_brts_density(pars1 = c(0.2,0.1,1,0.1), methode = "lsoda",brts = 1:10)
pbd_brts_density(pars1 = c(0.2,0.1,1,0.1), methode = "lsoda",brts = 1:10)

Cumulative density of duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_cumdensity computes the cumulative density of the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_cumdensity(pars, tau)
pbd_durspec_cumdensity(pars, tau)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[3]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`tau`	Value of the duration of speciation at which the cumulative density must be computed

Value

The cumulative density of the duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_cumdensity(pars = c(0.5,0.3,0.1),3)
pbd_durspec_cumdensity(pars = c(0.5,0.3,0.1),3)

Probability density for duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_density computes the probability density of the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_density(pars, tau)
pbd_durspec_density(pars, tau)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[3]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`tau`	The duration of speciation for which the density must be computed

Value

The probability density

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_density(pars = c(0.5,0.3,0.1), tau = 1)
pbd_durspec_density(pars = c(0.5,0.3,0.1), tau = 1)

Mean duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_mean computes the mean duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_mean(pars)
pbd_durspec_mean(pars)

Arguments

pars

Vector of parameters:

pars[1] corresponds to b (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate
pars[2] corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate
pars[3] corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species

Value

The expected duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_mean(pars = c(0.5,0.3,0.1)) 
pbd_durspec_mean(pars = c(0.5,0.3,0.1))

Actual calculation of the mean duration of speciation (equations 19 and 20 of reference article) assuming all inputs are correct

Description

Actual calculation of the mean duration of speciation (equations 19 and 20 of reference article) assuming all inputs are correct

Usage

pbd_durspec_mean_impl(la2, la3, mu2)
pbd_durspec_mean_impl(la2, la3, mu2)

Arguments

`la2`	lambda_2, the speciation completion rate, in probability per time unit
`la3`	lambda_3, speciation initiation rate of incipient species, in probability per time unit
`mu2`	mu_2 extinction rate of the incipient species, in probability per time unit

Author(s)

Rampal S. Etienne

References

Etienne, Rampal S., and James Rosindell. "Prolonging the past counteracts the pull of the present: protracted speciation can explain observed slowdowns in diversification." Systematic Biology 61.2 (2012): 204-213.

mode of the duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_mode computes the mode of the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_mode(pars)
pbd_durspec_mode(pars)

Arguments

pars

Value

The expected duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_mode(pars = c(0.5,0.3,0.1))
pbd_durspec_mode(pars = c(0.5,0.3,0.1))

Moments of duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_moment computes the moments of the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_moment(pars, order)
pbd_durspec_moment(pars, order)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[3]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`order`	order of the moment to compute (1 is first moment, giving the mean)

Value

The moment of the duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_moment(pars = c(0.5,0.3,0.1),2)
pbd_durspec_moment(pars = c(0.5,0.3,0.1),2)

Quantiles of duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_quantile computes a quantile of the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_quantile(pars, p)
pbd_durspec_quantile(pars, p)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[3]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`p`	Quantile (e.g. p = 0.5 gives the median)

Value

The quantil of the duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_quantile(pars = c(0.5,0.3,0.1),0.5)
pbd_durspec_quantile(pars = c(0.5,0.3,0.1),0.5)

Variance in duration of speciation under protracted birth-death model of diversification

Description

pbd_durspec_var computes the variance in the duration of speciation under the protracted speciation model for a given set of parameters

Usage

pbd_durspec_var(pars)
pbd_durspec_var(pars)

Arguments

pars

Value

The variance in the duration of speciation

Author(s)

Rampal S. Etienne

Examples

 pbd_durspec_var(pars = c(0.5,0.3,0.1))
pbd_durspec_var(pars = c(0.5,0.3,0.1))

Loglikelihood for protracted birth-death model of diversification

Description

pbd_loglik computes the loglikelihood of the parameters of the protracted speciation model given a set of branching times and number of missing species

Usage

pbd_loglik(
  pars1,
  pars1f = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  pars2 = c(1, 1, 2, 1, "lsoda", 0, 0),
  brts,
  missnumspec = 0
)
pbd_loglik(
  pars1,
  pars1f = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  pars2 = c(1, 1, 2, 1, "lsoda", 0, 0),
  brts,
  missnumspec = 0
)

Arguments

`pars1`	Vector of parameters: `pars1[1]` corresponds to b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `pars1[2]` corresponds to mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `pars1[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars1[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species When rates depend on time this time dependence should be specified in pars1f and pars1 then becomes the parameters used in pars1f
`pars1f`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `pars1f[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `pars1f[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `pars1f[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars1f[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`pars2`	Vector of model settings: `pars2[1]` set the conditioning on non-extinction of the clade (1) or not (0) `pars2[2]` sets whether the likelihood is for the branching times (0) or the phylogeny (1) `pars2[3]` sets whether the first element of the branching times is the stem (1) or the crown (2) age `pars2[4]` sets whether the parameters and likelihood should be shown on screen (1) or not (0) `pars2[5]` sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details. `pars2[6]`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) `pars2[7]`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)
`brts`	A set of branching times of a phylogeny, all positive
`missnumspec`	The number of species that are in the clade but missing in the phylogeny

Value

The loglikelihood

Author(s)

Rampal S. Etienne

Examples

 pbd_loglik(pars1 = c(0.2,0.1,1,0.1), pars2 = c(1,1,2,0,"lsoda"),brts = 1:10)
pbd_loglik(pars1 = c(0.2,0.1,1,0.1), pars2 = c(1,1,2,0,"lsoda"),brts = 1:10)

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-06, 1e-06, 1e-06),
  maxiter = 2000,
  optimmethod = "subplex",
  num_cycles = num_cycles,
  verbose = FALSE
)
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-06, 1e-06, 1e-06),
  maxiter = 2000,
  optimmethod = "subplex",
  num_cycles = num_cycles,
  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'.
`num_cycles`	Number of cycles of the optimization (default is 1).
`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` - peciation 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 econd 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

Calculate the mean duration of speciation (equations 19 and 20 of reference article), non-vectorized

Description

Calculate the mean duration of speciation (equations 19 and 20 of reference article), non-vectorized

Usage

pbd_mean_durspec(eri, scr, siri)
pbd_mean_durspec(eri, scr, siri)

Arguments

`eri`	one single extinction rate of the incipient species, or mu_2 in article, in probability per time unit
`scr`	one single speciation completion rate, or lambda_2 in article, in probability per time unit
`siri`	one single speciation initiation rate of incipient species, or lambda_3 in article, in probability per time unit

Value

the means duration of speciation, in time units

Author(s)

Richel J.C. Bilderbeek

References

Examples


  eri <- 0.1 # extinction rate of incipient species
  scr <- 0.2 # speciation completion rate
  siri <- 0.3 # speciation initiation rate of incipient species
  mean_durspec <- pbd_mean_durspec(eri, scr, siri)
  expected_mean_durspec <- 2.829762
  testthat::expect_equal(mean_durspec, expected_mean_durspec,
    tolerance = 0.000001)

eri <- 0.1 # extinction rate of incipient species
  scr <- 0.2 # speciation completion rate
  siri <- 0.3 # speciation initiation rate of incipient species
  mean_durspec <- pbd_mean_durspec(eri, scr, siri)
  expected_mean_durspec <- 2.829762
  testthat::expect_equal(mean_durspec, expected_mean_durspec,
    tolerance = 0.000001)

Calculate the mean durations of speciation (equations 19 and 20 of reference article)

Description

Calculate the mean durations of speciation (equations 19 and 20 of reference article)

Usage

pbd_mean_durspecs(eris, scrs, siris)
pbd_mean_durspecs(eris, scrs, siris)

Arguments

`eris`	one or more extinction rates of the incipient species, or mu_2 in article, in probability per time unit. These values will be recycled if needed
`scrs`	one or more speciation completion rates, or lambda_2 in article, in probability per time unit. These values will be recycled if needed
`siris`	one or more speciation initiation rates of incipient species, or lambda_3 in article, in probability per time unit. These values will be recycled if needed

Value

the means durations of speciation, in time units. Puts an NA at each invalid combination of inputs

Author(s)

Richel J.C. Bilderbeek

References

Examples


  eris <- c(0.1, 0.2) # extinction rates of incipient species
  scrs <- c(0.2, 0.3)  # speciation completion rates
  siris <- c(0.3, 0.4) # speciation initiation rates of incipient species
  mean_durspecs <- pbd_mean_durspecs(eris, scrs, siris)
  expected_mean_durspecs <- c(2.829762, 1.865386)
  testthat::expect_equal(mean_durspecs, expected_mean_durspecs,
    tolerance = 0.000001)

eris <- c(0.1, 0.2) # extinction rates of incipient species
  scrs <- c(0.2, 0.3)  # speciation completion rates
  siris <- c(0.3, 0.4) # speciation initiation rates of incipient species
  mean_durspecs <- pbd_mean_durspecs(eris, scrs, siris)
  expected_mean_durspecs <- c(2.829762, 1.865386)
  testthat::expect_equal(mean_durspecs, expected_mean_durspecs,
    tolerance = 0.000001)

Maximization of loglikelihood under protracted birth-death model of diversification

Description

Likelihood maximization for protracted birth-death model of diversification

Usage

pbd_ML(
  brts,
  initparsopt = c(0.2, 0.1, 1),
  idparsopt = 1:length(initparsopt),
  idparsfix = NULL,
  parsfix = NULL,
  exteq = 1,
  parsfunc = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  missnumspec = 0,
  cond = 1,
  btorph = 1,
  soc = 2,
  methode = "lsoda",
  n_low = 0,
  n_up = 0,
  tol = c(1e-06, 1e-06, 1e-06),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  optimmethod = "subplex",
  num_cycles = 1,
  verbose = TRUE
)
pbd_ML(
  brts,
  initparsopt = c(0.2, 0.1, 1),
  idparsopt = 1:length(initparsopt),
  idparsfix = NULL,
  parsfix = NULL,
  exteq = 1,
  parsfunc = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  missnumspec = 0,
  cond = 1,
  btorph = 1,
  soc = 2,
  methode = "lsoda",
  n_low = 0,
  n_up = 0,
  tol = c(1e-06, 1e-06, 1e-06),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  optimmethod = "subplex",
  num_cycles = 1,
  verbose = TRUE
)

Arguments

`brts`	A set of branching times of a phylogeny, all positive
`initparsopt`	The initial values of the parameters that must be optimized
`idparsopt`	The ids of the parameters that must be optimized, e.g. 1:4 for all parameters. The ids are defined as follows: id == 1 corresponds to b (speciation-initiation rate) id == 2 corresponds to mu_1 (extinction rate of good species) id == 3 corresponds to la_1 (speciation-completion rate) id == 4 corresponds to mu_2 (extinction rate of incipient species)
`idparsfix`	The ids of the parameters that should not be optimized, e.g. c(2,4) if mu_1 and mu_2 should not be optimized, but only b and la_1. In that case idparsopt must be c(1,3).
`parsfix`	The values of the parameters that should not be optimized
`exteq`	Sets whether incipient species have the same (1) or different (0) extinction rate as good species. If exteq = 0, then idparsfix and idparsopt should together have all parameters 1:4
`parsfunc`	Specifies functions how the rates depend on time, default functions are constant functions
`missnumspec`	The number of species that are in the clade but missing in the phylogeny
`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 number of extant taxa
`btorph`	Sets whether the likelihood is for the branching times (0) or the phylogeny (1)
`soc`	Sets whether the first element of the branching times is the stem (1) or the crown (2) age
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.
`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' (default of previous versions)
`num_cycles`	Number of cycles of the optimization (default is 1).
`verbose`	if TRUE, explanatory text will be shown

Value

A data frame with the following components:

`b`	gives the maximum likelihood estimate of b
`mu_1`	gives the maximum likelihood estimate of mu_1
`la_1`	gives the maximum likelihood estimate of la_1
`mu_2`	gives the maximum likelihood estimate of mu_2
`loglik`	gives the maximum loglikelihood
`df`	gives the number of estimated parameters, i.e. degrees of feedom
`conv`	gives a message on convergence of optimization; conv = 0 means convergence

Author(s)

Rampal S. Etienne

Examples


pbd_ML(1:10,initparsopt = c(4.640321,4.366528,0.030521), exteq = 1)

pbd_ML(1:10,initparsopt = c(4.640321,4.366528,0.030521), exteq = 1)

Mean number of species under protracted birth-death model of diversification

Description

pbd_numspec_mean computes the mean number of species under the protracted speciation model for a given set of parameters

Usage

pbd_numspec_mean(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda"
)
pbd_numspec_mean(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda"
)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_1 = la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to mu_1 (= mu_g in ER2012) = extinction rate of good species `pars[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`parsf`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `parsf[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `parsf[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `parsf[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `parsf[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`age`	the stem or crown age (see soc)
`soc`	specify whether it is the stem or the crown age
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.

Details

The number of representative species at time t = T when starting with one lineage at time t = 0 is geometrically distributed with parameter equal to g(T) where g(T) follows from the set of ODEs in Etienne et al. 2014.

Value

The expected number of representative species

Author(s)

Rampal S. Etienne

Examples

pbd_numspec_mean(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2)
pbd_numspec_mean(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2)

Median number of species under protracted birth-death model of diversification

Description

pbd_numspec_median computes the median number of species under the protracted speciation model for a given set of parameters

Usage

pbd_numspec_median(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda"
)
pbd_numspec_median(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda"
)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_1 = la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to mu_1 (= mu_g in ER2012) = extinction rate of good species `pars[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`parsf`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `parsf[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `parsf[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `parsf[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `parsf[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`age`	the stem or crown age (see soc)
`soc`	specify whether it is the stem or the crown age
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.

Value

The mean number of representative species

Author(s)

Rampal S. Etienne

Examples

pbd_numspec_median(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2)
pbd_numspec_median(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2)

Number of species for a given quantile under protracted birth-death model of diversification

Description

pbd_numspec_quantile computes the number of species for a given quantile under the protracted speciation model for a given set of parameters

Usage

pbd_numspec_quantile(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda",
  quantile
)
pbd_numspec_quantile(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  soc = 2,
  methode = "lsoda",
  quantile
)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_1 = la_3 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to mu_1 (= mu_g in ER2012) = extinction rate of good species `pars[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`parsf`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `parsf[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `parsf[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `parsf[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `parsf[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`age`	the stem or crown age (see soc)
`soc`	specify whether it is the stem or the crown age
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.
`quantile`	the quantile

Value

The number of species at a given quantile

Author(s)

Rampal S. Etienne

Examples

pbd_numspec_quantile(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2, quantile = 0.95)
pbd_numspec_quantile(pars = c(0.3,0.1,0.5,0.1), age = 10, soc = 2, quantile = 0.95)

Function to simulate the protracted speciation process

Description

Simulating the protracted speciation process using the Doob-Gillespie algorithm. This function differs from pbd_sim_cpp that 1) it does not require that the speciation-initiation rate is the same for good and incipient species, and 2) that it simulates the exact protracted speciation process, and not the approximation made by the coalescent point process. This function provides also the conversion to the approximation as output.

Usage

pbd_sim(
  pars,
  age,
  soc = 2,
  plotit = FALSE,
  limitsize = Inf,
  add_shortest_and_longest = FALSE
)
pbd_sim(
  pars,
  age,
  soc = 2,
  plotit = FALSE,
  limitsize = Inf,
  add_shortest_and_longest = FALSE
)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b_1, the speciation-initiation rate of good species `pars[2]` corresponds to la_1, the speciation-completion rate `pars[3]` corresponds to b_2, the speciation-initiation rate of incipient species `pars[4]` corresponds to mu_1, the extinction rate of good species `pars[5]` corresponds to mu_2, the extinction rate of incipient species
`age`	Sets the age for the simulation
`soc`	Sets whether this age is the stem (1) or crown (2) age
`plotit`	Sets whether the various trees produced by the function should be plotted or not
`limitsize`	Sets a maximum to the number of incipient + good species that are created during the simulation; if exceeded, the simulation is aborted and removed.
`add_shortest_and_longest`	Gives the output of the new samplemethods 'shortest' and 'longest'.

Value

out

A list with the following elements:

tree is the tree of extant species in phylo format
stree_random is a tree with one random sample per species in phylo format
stree_oldest is a tree with the oldest sample per species in phylo format
stree_youngest is a tree with the youngest sample per species in phylo format
L is a matrix of all events in the simulation where
- the first column is the incipient-level label of a species
- the second column is the incipient-level label of the parent of the species
- the third column is the time at which a species is born as incipient species
- the fourth column is the time of speciation-completion of the species
If the fourth element equals -1, then the species is still incipient. - the fifth column is the time of extinction of the species
If the fifth element equals -1, then the species is still extant. - The sixth column is the species-level label of the species
sL_random is a matrix like L but for stree_random
sL_oldest is a matrix like L but for stree_oldest
sL_youngest is a matrix like L but for stree_youngest
igtree.extinct is the tree in simmap format with incipient and good flags and including extinct species
igtree.extant is the tree in simmap format with incipient and good flags without extinct species
recontree is the reconstructed tree in phylo format, reconstructed using the approximation in Lambert et al. 2014
reconL is the matrix corresponding to recontree
L0 is a matrix where the crown age is at 0; for internal use only

Author(s)

Rampal S. Etienne

Examples

 pbd_sim(c(0.2,1,0.2,0.1,0.1),15)
pbd_sim(c(0.2,1,0.2,0.1,0.1),15)

Checked and safer version of PBD::pbd_sim

Description

Checked and safer version of PBD::pbd_sim

Usage

pbd_sim_checked(
  erg,
  eri,
  scr,
  sirg,
  siri,
  stem_age = NULL,
  crown_age = NULL,
  max_n_taxa = Inf,
  add_shortest_and_longest = FALSE
)
pbd_sim_checked(
  erg,
  eri,
  scr,
  sirg,
  siri,
  stem_age = NULL,
  crown_age = NULL,
  max_n_taxa = Inf,
  add_shortest_and_longest = FALSE
)

Arguments

`erg`	extinction rate of a good species
`eri`	extinction rate of an incipient species
`scr`	speciation completion rate
`sirg`	speciation initiation rate of a good species
`siri`	speciation initiation rate of an incipient species
`stem_age`	stem age. Set either the stem age or the crown age.
`crown_age`	crown age. Set either the stem age or the crown age.
`max_n_taxa`	maximum number of taxa. If this value is exceeded, the simulation is aborted and removed.
`add_shortest_and_longest`	Gives the output of the new samplemethods 'shortest' and 'longest'.

Function to simulate the approximate protracted speciation process

Description

Simulating the protracted speciation process according to the approxiomate model of Lambert et al. 2014. This function differs from pbd_sim that 1) it requires that the speciation-initiation rate is the same for good and incipient species, and 2) that it does not simulate the exact protracted speciation process, but an approximation made by the coalescent point process.

Usage

pbd_sim_cpp(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  ntips = NULL,
  soc = 2,
  plotltt = 1,
  methode = "lsoda"
)
pbd_sim_cpp(
  pars,
  parsf = c(function(t, pars) {     pars[1] }, function(t, pars) {     pars[2] },
    function(t, pars) {     pars[3] }, function(t, pars) {     pars[4] }),
  age,
  ntips = NULL,
  soc = 2,
  plotltt = 1,
  methode = "lsoda"
)

Arguments

`pars`	Vector of parameters: `pars[1]` corresponds to b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `pars[2]` corresponds to mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `pars[3]` corresponds to la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `pars[4]` corresponds to mu_2 (= mu_i in ER2012) = extinction rate of incipient species When rates depend on time this time dependence should be specified in pars1f and pars1 then becomes the parameters used in pars1f
`parsf`	Vector of functions how the rates depend on time, default functions are constant functions of the parameters in pars1: `parsf[1]` corresponds to time-dependence of b (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate `parsf[2]` corresponds to time-dependence of mu_1 (= mu_g in Etienne & Rosindell 2012) = extinction rate of good species `parsf[3]` corresponds to tiem-dependence of la_1 (= la_2 in Etienne & Rosindell 2012) = speciation completion rate `parsf[4]` corresponds to time-dependence of mu_2 (= mu_i in ER2012) = extinction rate of incipient species
`age`	Sets the crown age for the simulation
`ntips`	Set the number of tips. If NULL (default) the number of tips will be sampled
`soc`	Determines whether the simulation should start at stem (1) or crown (2) age
`plotltt`	Sets whether the lineage-through-time plot should be plotted (1) or not (0)
`methode`	Sets which method should be used in the ode-solver. Default is 'lsoda'. See package deSolve for details.

Value

A set of branching times

Author(s)

Rampal S. Etienne

Examples

 pbd_sim_cpp(pars = c(0.2,1,0.2,0.1),age = 15)
pbd_sim_cpp(pars = c(0.2,1,0.2,0.1),age = 15)

Package 'PBD'

Help Index

Protracted Birth-Death Model of Diversification

Description

Details

Author(s)

References

See Also

Bootstrap analysis under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Node depth probbaility density for protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Cumulative density of duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Probability density for duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Mean duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Actual calculation of the mean duration of speciation (equations 19 and 20 of reference article) assuming all inputs are correct

Description

Usage

Arguments

Author(s)

References

See Also

mode of the duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Moments of duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Quantiles of duration of speciation under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Variance in duration of speciation under protracted birth-death model of diversification

Description