Help for package PBD

Type:

Package

Title:

Protracted Birth-Death Model of Diversification

Version:

1.6

Date:

2026-09-04

Depends:

R (≥ 3.0.2)

Imports:

deSolve, ade4, ape, DDD, Matrix, phytools, graphics, stats, testit, utils

LinkingTo:

Rcpp (≥ 1.0.10), RcppEigen, BH (≥ 1.81.0-1)

Suggests:

knitr, rmarkdown, testthat

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>.

License:

GPL-2

Encoding:

UTF-8

VignetteBuilder:

knitr

URL:

https://github.com/rsetienne/PBD

BugReports:

https://github.com/rsetienne/PBD/issues

RoxygenNote:

7.3.3

NeedsCompilation:

yes

Packaged:

2026-09-04 15:59:44 UTC; rampa

Author:

Rampal S. Etienne

[aut, cre], Bart Haegeman

[aut], Richèl J.C. Bilderbeek

[ctb], Thijs Janzen

[ctb], Raphael Scherrer

[ctb]

Maintainer:

Rampal S. Etienne <r.s.etienne@rug.nl>

Repository:

CRAN

Date/Publication:

2026-09-04 16:40:02 UTC

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
)

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

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
)

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

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

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)

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
)

Arguments

brts

A set of branching times of a phylogeny, all positive

initparsopt

The initial values of the parameters that must be optimized

idparsopt

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

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
)

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)

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)

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)

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)

Arguments

pars

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)

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)

Arguments

pars

Value

The expected duration of speciation

Author(s)

Rampal S. Etienne

Examples

 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)

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)

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))

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)

Arguments

pars

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)

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)

Arguments

pars

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)

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)

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))

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
)

Arguments

pars1

pars1f

pars2

Vector of model settings:

pars2[1] sets 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)

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)

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)

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)

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)

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"
)

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)

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"
)

Arguments

pars

parsf

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)

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
)

Arguments

pars

parsf

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)

Mean number of good and incipient species under protracted birth-death model of diversification

Description

pbd_numspec_vec computes the mean number of good and incipient species under the protracted speciation model for a given set of parameters

Usage

pbd_numspec_vec(pars, age, initvec = c(1, 0))

Arguments

pars

Vector of parameters:

pars[1] corresponds to b_G (= la_1 in Etienne & Rosindell R2012) = speciation initiation rate of good species
pars[2] corresponds to b_I (= la_3 in Etienne & Rosindell R2012) = speciation initiation rate of incipient species
pars[3] corresponds to mu_G = extinction rate of good species
pars[4] corresponds to mu_I = extinction rate of incipient species
pars[5] corresponds to la (= la_2 in Etienne & Rosindell 2012) = speciation completion rate

age

the stem age

initvec

initial vector of good and incipient species. Default is 1 good species.

Details

The result follows from a simple 2D system of differential equations for G and I.

Value

The expected number of good and incipient species

Author(s)

Rampal S. Etienne

Examples

pbd_numspec_vec(pars = c(0.3,0.3,0.1,0.1,0.2), age = 10)

Mean number of good species with n incipient species and of n orphaned incipient species under protracted birth-death model of diversification

Description

pbd_numspec_vec2 computes the mean number of good species with n incipient daughter species and the number of cases with n orphaned incipient siblings under the protracted speciation model for a given set of parameters

Usage

pbd_numspec_vec2(
  pars = c(0.3, 0.3, 0.1, 0.1, 100),
  age = 10,
  initvec = c(1, rep(0, 20)),
  abstolint = 1e-10,
  reltolint = 1e-10,
  methode = "odeint::runge_kutta_cash_karp54"
)

Arguments

pars

age

the stem age

initvec

initial vector. Default is 1 good species with 0 incipient species.

abstolint

absolute tolerance in integration

reltolint

relative tolerance in integration

methode

integration method, one of the odeint solvers

Value

The expected number of good species with n incipient specis and the number of caes with n orphaned incipient siblings

Author(s)

Rampal S. Etienne

Examples

pbd_numspec_vec2(pars = c(0.3,0.3,0.1,0.1,0.2), age = 10)

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
)

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)

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
)

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"
)

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

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)

Package {PBD}

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

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Maximization of loglikelihood under protracted birth-death model of diversification

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

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