getBestShiftConfiguration {BAMMtools}

Get the best (sampled) rate shift configuration from a BAMM analysis

Description

Get the rate shift configuration with the maximum a posteriori probability, e.g., the shift configuration that was sampled most frequently with BAMM.

Usage

getBestShiftConfiguration(x, expectedNumberOfShifts, threshold = 5)

Arguments

x

Either a bammdata object or a credibleshiftset object.

expectedNumberOfShifts

The expected number of shifts under the prior.

threshold

The marginal posterior-to-prior odds ratio used as a cutoff for distinguishing between "core" and "non-core" rate shifts.

Details

This function estimates the rate shift configuration with the highest maximum a posteriori (MAP) probability. It returns a bammdata object with a single sample. This can be plotted with plot.bammdata, and individual rate shifts can then be added with addBAMMshifts.

The parameters of this object are averaged over all samples in the posterior that were assignable to the MAP shift configuration. All non-core shifts have been excluded, such that the only shift information contained in the object is from the "significant" rate shifts, as determined by the relevant marginal posterior-to-prior odds ratio threshold.

You can extract the same information from the credible set of shift configurations. See credibleShiftSet for more information.

Value

A class bammdata object with a single sample, corresponding to the diversification rate shift configuration with the maximum a posteriori probability. See getEventData for details.

Author(s)

Dan Rabosky

See Also

getEventData, credibleShiftSet, plot.credibleshiftset, plot.bammdata

Examples

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

# Get prior distribution on shifts-per-branch:
bp <- getBranchShiftPriors(whales, expectedNumberOfShifts = 1)

# Pass the event data object in to the function:
best <- getBestShiftConfiguration(ed, expectedNumberOfShifts = 1,
                                  threshold = 5)
plot(best, lwd=2)
addBAMMshifts(best, cex=2)

# Now we can also work with the credible shift set:
css <- credibleShiftSet(ed, expectedNumberOfShifts = 1, threshold = 5)

summary(css)

# examine model-averaged shifts from MAP configuration-
# This gives us parameters, times, and associated nodes
#   of each evolutionary rate regime (note that one of
#   them corresponds to the root)
css$eventData[[1]];

# Get bammdata representation of MAP configuration:
best <- getBestShiftConfiguration(css, expectedNumberOfShifts = 1,
                                  threshold = 5)

plot(best)
addBAMMshifts(best)