Calculates the likelihood of a phylogeny exactly as is done
by BAMM, given a set of events.
Usage
BAMMlikelihood(phy, eventdata, gen = "last", segLength = 0.02, sf = 1,
return.intermediates = FALSE, e_prob_condition = "if_different", ...)
Arguments
phy
Either an object of class phylo or the path to a tree
file in newick format.
eventdata
A table of event data, as returned by BAMM, either
as an object of class dataframe or as the path to an event_data
file. Alternatively, a named numeric vector of length two holding
speciation ("lambda") and extinction ("mu") rates for the
constant-rate birth-death process.
gen
The BAMM generation for which the likelihood should be
calculated. Can be an integer specifying a specific generation, or
last, specifying the last generation, or all, in which
case the likelihood will be calculated for all generations.
segLength
The relative segment length, exactly as defined for
BAMM.
sf
The sampling fraction.
return.intermediates
Debugging option, returns augmented
phylo objects for each generation, see Details.
e_prob_condition
Approach for how extinction probabilities are
handled at nodes.
...
Additional arguments that will be passed to an internal
function computeBAMMlikelihood.
Details
This function allows the user to check the likelihoods computed
by BAMM using an independent R-based implementation. This is
designed to provide a check on potential software bugs that might be
introduced during future BAMM development and which might
compromise the likelihood calculation. If you observe measurable
discrepancies between the likelihood computed by this function and the
corresponding likelihood returned by BAMM, please inform the
BAMM development team.
Value
If return.intermediates == TRUE, then phylo objects
are returned with the following components:
event_times
A list of length (number of nodes), where
event_times[[k]] is the vector of absolute times, in order, of
events that happened on a focal branch. If no event, it is
NULL.
event_id
A list of length equal to number of nodes, as
event_times, but holding the corresponding event id.
events
A dataframe giving parameters and associated nodes (and
unique index values) of the event data.
node_event
The event governing the process realized at the
node. This will be the first event encountered as one moves
rootwards towards the tips from the focal node.
Author(s)
Dan Rabosky, Pascal Title
Examples
# a global sampling fraction of 0.98 was used in generating the whales
# dataset.
data(whales, events.whales, mcmc.whales)
x <- BAMMlikelihood(whales, events.whales, gen = 'last', sf = 0.98)
# Does the likelihood generated by BAMM match the R implementation?
identical(round(x, 3), mcmc.whales[nrow(mcmc.whales), 'logLik'])
# an example with a constant-rate birth-death process:
pars <- c(0.5, 0.45)
names(pars) <- c("lambda", "mu")
BAMMlikelihood(whales, pars, sf = 0.98)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(BAMMtools)
Loading required package: ape
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/BAMMtools/BAMMlikelihood.Rd_%03d_medium.png", width=480, height=480)
> ### Name: BAMMlikelihood
> ### Title: Calculate 'BAMM' likelihood
> ### Aliases: BAMMlikelihood
>
> ### ** Examples
>
> # a global sampling fraction of 0.98 was used in generating the whales
> # dataset.
> data(whales, events.whales, mcmc.whales)
>
> x <- BAMMlikelihood(whales, events.whales, gen = 'last', sf = 0.98)
>
> # Does the likelihood generated by BAMM match the R implementation?
> identical(round(x, 3), mcmc.whales[nrow(mcmc.whales), 'logLik'])
[1] TRUE
>
> # an example with a constant-rate birth-death process:
> pars <- c(0.5, 0.45)
> names(pars) <- c("lambda", "mu")
> BAMMlikelihood(whales, pars, sf = 0.98)
[1] -311.4387
>
>
>
>
>
> dev.off()
null device
1
>