loglikelihood_hbd {castor} | R Documentation |

Given a rooted ultrametric timetree, and a homogenous birth-death (HBD) model, i.e., with speciation rate *λ*, extinction rate *μ* and sampling fraction *ρ*, calculate the likelihood of the tree under the model. The speciation and extinction rates may be time-dependent. “Homogenous” refers to the assumption that, at any given moment in time, all lineages exhibit the same speciation/extinction rates (in the literature this is sometimes referred to simply as “birth-death model”). Alternatively to *λ* and *μ*, the likelihood may also be calculated based on the pulled diversification rate (PDR; Louca et al. 2018) and the product *ρ(0)\cdotλ(0)*, or based on the pulled speciation rate (PSR). In either case, the time-profiles of *λ*, *μ*, the PDR or the PSR are specified as piecewise polynomially functions (splines), defined on a discrete grid of ages.

loglikelihood_hbd(tree, oldest_age = NULL, age0 = 0, rho0 = NULL, rholambda0 = NULL, age_grid = NULL, lambda = NULL, mu = NULL, PDR = NULL, PSR = NULL, splines_degree = 1, condition = "auto", max_model_runtime = -1, relative_dt = 1e-3)

`tree` |
A rooted ultrametric tree of class "phylo". |

`oldest_age` |
Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then |

`age0` |
Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which |

`rho0` |
Numeric between 0 (exclusive) and 1 (inclusive), specifying the sampling fraction of the tree at |

`rholambda0` |
Strictly positive numeric, specifying the product of the sampling fraction and the speciation rateat |

`age_grid` |
Numeric vector, listing discrete ages (time before present) on which either |

`lambda` |
Numeric vector, of the same size as |

`mu` |
Numeric vector, of the same size as |

`PDR` |
Numeric vector, of the same size as |

`PSR` |
Numeric vector, of the same size as |

`splines_degree` |
Integer, either 0,1,2 or 3, specifying the polynomial degree of the provided |

`condition` |
Character, either "crown", "stem", "auto" or "none" (the last one is only available if |

`max_model_runtime` |
Numeric, maximum allowed runtime (in seconds) for evaluating the likelihood. If the likelihood calculation takes longer than this (appoximate) threshold, it halts and returns with an error. If negative (default), this option is ignored. |

`relative_dt` |
Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient. |

If `age0>0`

, the input tree is essentially trimmed at `age0`

(omitting anything younger than `age0`

), and the is likelihood calculated for the trimmed tree while shifting time appropriately. In that case, `rho0`

is interpreted as the sampling fraction at `age0`

, i.e. the fraction of lineages extant at `age0`

that are repreented in the tree. Similarly, `rholambda0`

is the product of the sampling fraction and *λ* at `age0`

.

This function supports three alternative parameterizations of HBD models, either using the speciation and extinction rates and sampling fraction (*λ*, *μ* and *ρ(τ_o)* (for some arbitrary age *τ_o*), or using the pulled diversification rate (PDR) and the product *ρ(τ_o)\cdotλ(τ_o* (sampling fraction times speciation rate at *τ_o*), or using the pulled speciation rate (PSR). The latter two options should be interpreted as a parameterization of congruence classes, i.e. sets of models that have the same likelihood, rather than specific models, since multiple combinations of *λ*, *μ* and *ρ(τ_o)* can have identical PDRs, *ρ(τ_o)\cdotλ(τ_o)* and PSRs (Louca and Pennell, in review).

For large trees the asymptotic time complexity of this function is O(Nips). The tree may include monofurcations as well as multifurcations, and the likelihood formula accounts for those (i.e., as if monofurcations were omitted and multifurcations were expanded into bifurcations).

A named list with the following elements:

`success` |
Logical, indicating whether the calculation was successful. If |

`loglikelihood` |
Numeric. If |

Stilianos Louca

H. Morlon, T. L. Parsons, J. B. Plotkin (2011). Reconciling molecular phylogenies with the fossil record. Proceedings of the National Academy of Sciences. 108:16327-16332.

S. Louca et al. (2018). Bacterial diversification through geological time. Nature Ecology & Evolution. 2:1458-1467.

S. Louca and M. W. Pennell (in review as of 2019)

# generate a random tree with constant rates Ntips = 100 params = list(birth_rate_factor=1, death_rate_factor=0.2, rarefaction=0.5) tree = generate_random_tree(params, max_tips=Ntips, coalescent=TRUE)$tree # get the loglikelihood for an HBD model with the same parameters that generated the tree # in particular, assuming time-independent speciation & extinction rates LL = loglikelihood_hbd( tree, rho0 = params$rarefaction, age_grid = NULL, # assume time-independent rates lambda = params$birth_rate_factor, mu = params$death_rate_factor) if(LL$success){ cat(sprintf("Loglikelihood for constant-rates model = %g\n",LL$loglikelihood)) } # get the likelihood for a model with exponentially decreasing (in forward time) lambda & mu beta = 0.01 # exponential decay rate of lambda over time age_grid = seq(from=0, to=100, by=0.1) # choose a sufficiently fine age grid lambda = 1*exp(beta*age_grid) # define lambda on the age grid mu = 0.2*lambda # assume similarly shaped but smaller mu LL = loglikelihood_hbd( tree, rho0 = params$rarefaction, age_grid = age_grid, lambda = lambda, mu = mu) if(LL$success){ cat(sprintf("Loglikelihood for exponential-rates model = %g\n",LL$loglikelihood)) }

[Package *castor* version 1.6.8 Index]