credibleShiftSet: Credible set of macroevolutionary rate shift configurations from BAMM results

Description

Computes the 95% (or any other %) credible set of macroevolutionary rate shift configurations from a bammdata object. These results can be analyzed further and/or plotted.

Usage

credibleShiftSet(ephy, expectedNumberOfShifts, threshold = 5, set.limit = 0.95, ...)

Arguments

ephy

an object of class bammdata

expectedNumberOfShifts

The expected number of rate shifts under the prior

threshold

The marginal posterior-to-prior odds ratio for a rate shift on a specific branch, used to distinguish core and non-core shifts

set.limit

the desired limit to the credible set. A value of 0.95 will return the 95% credible set of shift configurations.

...

Other arguments to credibleShiftSet

Value

A class credibleshiftset object with many components. Most components are an ordered list of length L, where L is the number of distinct shift configurations in the credible set. The first list element in each case corresponds to the shift configuration with the maximum a posteriori probability.
frequencyA vector of frequencies of shift configurations, including those that account for set.limit (typically, 0.95 or 0.99) of the probability of the data. The index of the i'th element of this vector is the i'th most probable shift configuration (excepting ties).
shiftnodesA list of the "core" rate shifts (marginal probability > threshold) that occurred in each distinct shift configuration in the credible set. The i'th vector from this list gives the core shift nodes for the i'th shift configuration. They are sorted by frequency, so x$shiftnodes[[1]] gives the shift nodes that occurred together in the shift configuration with the highest posterior probability.
indicesA list of vectors containing the indices of samples in the bammdata object that are assigned to a given shift configuration. All are sorted by frequency.
cumulativeLike frequency, but contains the cumulative frequencies
thresholdThe marginal posterior-to-prior odds for rate shifts on branches used during enumeration of distinct shift configurations.
number.distinctNumber of distinct shift configurations in the credible set
set.limitwhich credible set is this (0.9, 0.95, etc)?
coreshiftsA vector of node numbers corresponding to the core shifts. All of these nodes have a Bayes factor of at least BFcriterion supporting a rate shift.
In addition, a number of components that are defined similarly in class phylo or class bammdata objects:
edgeSee documentation for class phylo in package ape.
NnodeSee documentation for class phylo in package ape.
tip.labelSee documentation for class phylo in package ape.
edge.lengthSee documentation for class phylo in package ape.
beginThe beginning time of each branch in absolute time (the root is set to time zero)
endThe ending time of each branch in absolute time.
numberEventsAn integer vector with the number of core events contained in the bammdata object for each shift configuration in the credible set. The length of this vector is equal to the number of distinct shift configurations in the credible set.
eventDataA list of dataframes. Each element holds the average rate and location parameters for all samples from the posterior that were assigned to a particular distinct shift configuration. Each row in a dataframe holds the data for a single event. Data associated with an event are: node - a node number. This identifies the branch where the event originates. time - this is the absolute time on that branch where the event originates (with the root at time 0). lam1 - an initial rate of speciation or trait evolution. lam2 - a decay/growth parameter. mu1 - an initial rate of extinction. mu2 - a decay/growth parameter. index - a unique integer associated with the event. See 'Details' in the documentation for getEventData for more information.
eventVectorsA list of integer vectors. Each element is for a single shift configuration in the posterior. For each branch in the bammdata object, gives the index of the event governing the (tipwards) end of the branch. Branches are ordered increasing here and elsewhere.
eventBranchSegsA list of matrices. Each element of the list is a single distinct shift configuration. Each matrix has four columns: Column 1 identifies a node in phy. Column 2 identifies the beginning time of the branch or segment of the branch that subtends the node in Column 1. Column 3 identifies the ending time of the branch or segment of the branch that subtends the node in Column 1. Column 4 identifies the index of the event that occurs along the branch or segment of the branch that subtends the node in Column 1.
tipStatesA list of integer vectors. Each element is a single distinct shift configuration. For each tip the index of the event that occurs along the branch subtending the tip. Tips are ordered increasing here and elsewhere.
tipLambdaA list of numeric vectors. Each element is a single distinct shift configuration. For each tip is the average rate of speciation or trait evolution at the end of the terminal branch subtending that tip (averaged over all samples that are assignable to this particular distinct shift configuration).
tipMuA list of numeric vectors. Each element is a single distinct shift configuration. For each tip the rate of extinction at the end of the terminal branch subtending that tip. Meaningless if working with BAMM trait results.
typeA character string. Either "diversification" or "trait" depending on your BAMM analysis.
downseqAn integer vector holding the nodes of phy. The order corresponds to the order in which nodes are visited by a pre-order tree traversal.
lastvisitAn integer vector giving the index of the last node visited by the node in the corresponding position in downseq. downseq and lastvisit can be used to quickly retrieve the descendants of any node. e.g. the descendants of node 89 can be found by downseq[which(downseq==89):which(downseq==lastvisit[89]).

Details

Computes the 95% credible set (or XX% credible set, depending on set.limit) of diversification shift configurations sampled using BAMM. This is analogous to a credible set of phylogenetic tree topologies from a Bayesian phylogenetic analysis.

To understand how this calculation is performed, one must first distinguish between "core" and "non-core" rate shifts. A "core shift" is a rate shift with a marginal probability that is substantially elevated above the probability expected on the basis of the prior alone. With BAMM, every branch in a phylogenetic tree is associated with some non-zero prior probability of a rate shift. Typically this is a very low per-branch shift probability (this prior is determined by the value of the "poissonRatePrior" parameter in a BAMM analysis).

If we compute distinct shift configurations with every sampled shift (including those shifts with very low marginal probabilities), the number of distinct shift configurations will be overwhelmingly high. However, most of these configurations include shifts with marginal probabilities that are expected even under the prior alone. Hence, using these shifts to identify distinct shift configurations simply generates noise and isn't particularly useful.

The solution adopted in BAMMtools is, for each branch in the phylogeny, to compute both the posterior and prior probabilities of a rate shift occurring. The ratio of these probabilities is a branch-specific marginal odds ratio: it is the marginal posterior frequency of one or more rate shifts normalized by the corresponding prior probability. Hence, any branch with a marginal odds ratio of 1.0 is one where the observed (posterior) odds of a rate shift are no different from the prior odds. A value of 10 implies that the posterior probability is 10 times the prior probability.

The user of credibleShiftSet must specify a threshold argument. This is simply a cutoff value for identifying "important" shifts for the purposes of identifying distinct shift configurations. This does not imply that it is identifying "significant" shifts. See the online documentation on this topic available at www.bamm-project.org for more information. If you specify threshold = 5 as an argument to credibleShiftSet, the function will ignore all branches with marginal odds ratios less than 5 during the enumeration of topologically distinct shift configurations. Only shifts with marginal odds ratios greater than or equal to threshold will be treated as core shifts for the purposes of identifying distinct shift configurations.

For each shift configuration in the credible set, this function will compute the average diversification parameters. For example, the most frequent shift configuration (the maximum a posteriori shift configuration) might have 3 shifts, and 150 samples from your posterior (within the bammdata object) might show this shift configuration. However, the parameters associated with each of these shift configurations (the actual evolutionary rate parameters) might be different for every sample. This function returns the mean set of rate parameters for each shift configuration, averaging over all samples from the posterior that can be assigned to a particular shift configuration.

See the examples below for a worked example.

References

bamm-project.org

Examples

Run this code

data(events.whales, whales)
ed <- getEventData(whales, events.whales, burnin=0.1, nsamples=500)

cset <- credibleShiftSet(ed, expectedNumberOfShifts = 1, threshold = 5)

# Here is the total number of samples in the posterior:
length(ed$eventData)

# And here is the number of distinct shift configurations:
cset$number.distinct

# here is the summary statistics:
summary(cset)

# Accessing the raw frequency vector for the credible set:
cset$frequency

#The cumulative frequencies:
cset$cumulative

# The first element is the shift configuration with the maximum 
#  a posteriori probability. We can identify all the samples from 
#  posterior that show this shift configuration:

cset$indices[[1]]
 
# Now we can plot the credible set:
plot.credibleshiftset(cset, plotmax=4)

Run the code above in your browser using DataLab