Package 'RevGadgets'

Description

This function computes the Bayes Factor in favor of a rate-shift between time t1 and t2 (t1 < t2). The default assumption (suitable to standard HSMRF and GMRF models) is that the prior probability of a shift is 0.5.

Usage

calculateShiftBayesFactor(
  rate_trace,
  time_trace,
  rate_name,
  time_name,
  t1,
  t2,
  prior_prob = 0.5,
  decrease = TRUE,
  return_2lnBF = TRUE
)

Arguments

Value

The Bayes Factor.

References

Kass and Raftery (1995) Bayes Factors. JASA, 90 (430), 773-795.

Examples

#' # download the example datasets to working directory
url_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_times.log"
dest_path_times <- "primates_EBD_speciation_times.log"
download.file(url_times, dest_path_times)

url_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_rates.log"
dest_path_rates <- "primates_EBD_speciation_rates.log"
download.file(url_rates, dest_path_rates)

# to run on your own data, change this to the path to your data file
speciation_time_file <- dest_path_times
speciation_rate_file <- dest_path_rates

speciation_times <- readTrace(speciation_time_file, burnin = 0.25)
speciation_rate <- readTrace(speciation_rate_file, burnin = 0.25)

calculateShiftBayesFactor(speciation_rate,
                          speciation_times,
                          "speciation",
                          "interval_times",
                          0.0,40.0,
                          decrease=FALSE)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_times, dest_path_rates)

Description

Produce default RevGadgets colors

Usage

colFun(n)

Arguments

Details

Produces a vector of colors from the default RevGadgets colors of length given by n, maximum of 12 colors.

Value

Character vector of color hex codes.

Examples

my_colors <- colFun(2)

Description

Combine traces into one trace file

Usage

combineTraces(traces, burnin = 0)

Arguments

Details

Combines multiple traces from independent MCMC replicates into one trace file.

Value

combineTraces() returns a list of data frames of length 1, corresponding to the combination of the provided samples.

Examples

#' # download the example dataset to working directory
url_1 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_1.log"
dest_path_1 <- "primates_cytb_GTR_run_1.log"
download.file(url_1, dest_path_1)

url_2 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_2.log"
dest_path_2 <- "primates_cytb_GTR_run_2.log"
download.file(url_2, dest_path_2)

# to run on your own data, change this to the path to your data file
file_1 <- dest_path_1
file_2 <- dest_path_2

# read in the multiple trace files
multi_trace <- readTrace(path = c(file_1, file_2), burnin = 0.0)

# combine samples after discarding 10% burnin
combined_trace <- combineTraces(trace = multi_trace,
                                burnin = 0.1)

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_1, dest_path_2)

Description

This function plots a distribution of trees (e.g obtained from an MCMC inference) with branch-specific rates or other data. The plot is similar to those produced by DensiTree, i.e all the trees are overlapped with each other. The data is expected to be given per node, and will be associated with the branch above its corresponding node. Its values are plotted as a color gradient.

Usage

densiTreeWithBranchData(
  tree_files = NULL,
  burnin = 0.1,
  trees = NULL,
  data = NULL,
  data_name = NULL,
  type = "cladogram",
  consensus = NULL,
  direction = "rightwards",
  scaleX = FALSE,
  width = 1,
  lty = 1,
  cex = 0.8,
  font = 3,
  tip.color = 1,
  adj = 0,
  srt = 0,
  keep_underscores = FALSE,
  label_offset = 0.01,
  scale_bar = TRUE,
  jitter = list(amount = 0, random = TRUE),
  color_gradient = c("red", "yellow", "green"),
  alpha = NULL,
  bias = 1,
  data_intervals = NULL,
  ...
)

Arguments

Details

If no consensus tree is provided, a consensus tree will be computed. This should avoid too many unnecessary crossings of edges. Trees should be rooted, other wise the output may not be visually pleasing. The jitter parameter controls whether to shift trees so that they are not exactly on top of each other. If amount = 0, no jitter is applied. If random = TRUE, the applied jitter is calculated as runif(n, -amount, amount), otherwise seq(-amount, amount, length=n), where n is the number of trees.

Value

No return value, produces plot in base R

References

This code is adapted from the densiTree function by Klaus Schliep [email protected]. densiTree is inspired from the DensiTree program by Remco Bouckaert.

Remco R. Bouckaert (2010) DensiTree: making sense of sets of phylogenetic trees Bioinformatics, 26 (10), 1372-1373.

Examples

# generate random trees & data
trees <- lapply(1:5, function(x) ape::rcoal(5))
data <- lapply(1:5, function(x) stats::runif(9, 1, 10))

# densiTree plot
densiTreeWithBranchData(trees = trees, data = data, width = 2)

# densiTree plot with different colors
densiTreeWithBranchData(trees = trees, data = data,
                        color_gradient = c("green", "blue"), width = 2)

Description

Drop one or multiple tips from your tree

Usage

dropTip(tree, tips)

Arguments

Details

Modifies a tree object (in RevGadget's format) by dropping one or more tips from the tree and from any associated data. Wrapper for treeio::drop.tip().

Value

returns a list of list of treedata objects, with the modified tips.

See Also

treeio: drop.tip and ape: drop.tip.

Examples

file <- system.file("extdata",
                    "sub_models/primates_cytb_GTR_MAP.tre",
                    package="RevGadgets")
tree <- readTrees(paths = file)
tree_dropped <- dropTip(tree, "Otolemur_crassicaudatus")

Description

Modified from RmcdrPlugin.KMggplot2 step ribbon plots.

Usage

geom_stepribbon(
  mapping = NULL,
  data = NULL,
  stat = "identity",
  position = "identity",
  na.rm = FALSE,
  show.legend = NA,
  inherit.aes = TRUE,
  ...
)

Arguments

Details

geom_stepribbon is an extension of the geom_ribbon, and is optimized for Kaplan-Meier plots with pointwise confidence intervals or a confidence band.

See Also

geom_ribbon geom_stepribbon inherits from geom_ribbon. geom_stepribbon is modified from RcmdrPlugin.KMggplot2::geom_stepribbon.

Examples

huron <- data.frame(year = 1875:1972, level = as.vector(LakeHuron))

h <- ggplot2::ggplot(huron, ggplot2::aes(year))

h + geom_stepribbon(ggplot2::aes(ymin = level - 1, ymax = level + 1),
                    fill = "grey70") +
    ggplot2::geom_step(ggplot2::aes(y = level))

# contrast ggplot2::geom_ribbon with geom_stepribbon:
h + ggplot2::geom_ribbon(ggplot2::aes(ymin = level - 1, ymax = level + 1),
                         fill = "grey70") +
    ggplot2::geom_line(ggplot2::aes(y = level))

Description

Calculates the Maximum a Posteriori estimate for the trace of a quantitative variable

Usage

getMAP(var)

Arguments

Details

Uses the SANN method of the optim() function to approximate the MAP estimate

Value

the MAP estimate

See Also

optim

Examples

# download the example dataset to working directory
url <-
  "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.log"
dest_path <- "primates_cytb_GTR.log"
download.file(url, dest_path)

# to run on your own data, change this to the path to your data file
file <- dest_path

trace <- readTrace(paths = file)
MAP <- getMAP(trace[[1]]$"pi[1]")

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path)

Description

match Nodes

Usage

matchNodes(phy)

Arguments

Value

a data frame that translates ape node numbers to RevBayes node numbers

Examples

treefile <- system.file("extdata", "bds/primates.tre", package="RevGadgets")
tree <- readTrees(treefile)
map <- matchNodes(tree[[1]][[1]]@phylo)

Description

Plots the MAP estimates of ancestral states. Can accommodate cladogenetic reconstructions by plotting on shoulders. Defaults to varying the symbols by color to indicate estimated ancestral state and varying the size of the symbol to indicate the posterior probability of that estimate, but symbol shape may also vary to accommodate black and white figures. For more details on the aesthetics options, see parameter details below. For data with many character states (such as chromosome counts), vary the size of the symbol by estimated ancestral state, and vary the posterior probability of that estimate by a color gradient. Text labels at nodes and tips are also available.

Usage

plotAncStatesMAP(
  t,
  cladogenetic = FALSE,
  tip_labels = TRUE,
  tip_labels_size = 2,
  tip_labels_offset = 1,
  tip_labels_italics = FALSE,
  tip_labels_formatted = FALSE,
  tip_labels_remove_underscore = TRUE,
  tip_labels_states = FALSE,
  tip_labels_states_size = 2,
  tip_labels_states_offset = 0.1,
  node_labels_as = NULL,
  node_labels_size = 2,
  node_labels_offset = 0.1,
  node_labels_centered = FALSE,
  node_size_as = "state_posterior",
  node_color_as = "state",
  node_shape_as = NULL,
  node_shape = 19,
  node_color = "default",
  node_size = c(2, 6),
  tip_states = TRUE,
  tip_states_size = node_size,
  tip_states_shape = node_shape,
  state_transparency = 0.75,
  tree_layout = "rectangular",
  timeline = FALSE,
  geo = timeline,
  geo_units = list("epochs", "periods"),
  time_bars = timeline,
  ...
)

Arguments

Value

A ggplot object

Examples

# Standard ancestral state reconstruction example with various aesthetics

# process file
file <- system.file("extdata",
                    "comp_method_disc/ase_freeK.tree",
                    package="RevGadgets")
example <- processAncStates(file,
                            state_labels = c("1" = "Awesome",
                                             "2" = "Beautiful",
                                             "3" = "Cool!"))

# have states vary by color and indicate state pp with size (default)
plotAncStatesMAP(t = example)

# have states vary by color and indicate state pp with size ,
# and add a timeline
plotAncStatesMAP(t = example, timeline = TRUE)

# have states vary by color and symbol, label nodes with pp of states
plotAncStatesMAP(t = example,  node_shape_as = "state",
                 node_size = 4, node_shape = c(15, 17,20),
                 node_size_as = NULL, node_labels_as = "state_posterior")

# black and white figure - state as symbols and state pp with text
plotAncStatesMAP(t = example, node_color_as = NULL,
                 node_shape_as = "state", node_shape =  c(15, 17,20),
                 node_size_as = NULL, node_size = 4,
                 node_labels_as = "state_posterior",
                 node_color = "grey", state_transparency = 1)

# default with circular tree
plotAncStatesMAP(t = example, tree_layout = "circular")


# Chromosome evolution example

# process file
file <- system.file("extdata",
                    "chromo/ChromEvol_simple_final.tree",
                    package="RevGadgets")
chromo_example <- processAncStates(file, labels_as_numbers = TRUE)

# plot
plotAncStatesMAP(t = chromo_example, node_color_as = "state_posterior",
                 node_size_as = "state", node_color = colFun(2),
                 tip_labels_offset = 0.005, node_labels_as = "state",
                 node_labels_offset = 0, tip_labels_states = TRUE,
                 tip_labels_states_offset = 0, tip_states = FALSE)

# DEC example (cladogenetic)

# process file
file <- system.file("extdata", "dec/simple.ase.tre", package="RevGadgets")
labs <- c("1" = "K", "2" = "O", "3" = "M", "4" = "H", "5" = "KO",
"6" = "KM", "7" = "OM", "8" = "KH", "9" = "OH", "10" = "MH", "11" = "KOM",
"12" = "KOH", "13" = "KMH", "14" = "OMH", "15" = "KOMH")
dec_example <- processAncStates(file, state_labels = labs)

# plot
plotAncStatesMAP(t = dec_example,
                 cladogenetic = TRUE,
                 tip_labels_offset = 0.5)

Description

Plot character states and posterior probabilities as pies on nodes.

Usage

plotAncStatesPie(
  t,
  cladogenetic = FALSE,
  tip_labels = TRUE,
  tip_labels_size = 2,
  tip_labels_offset = 1,
  tip_labels_italics = FALSE,
  tip_labels_formatted = FALSE,
  tip_labels_remove_underscore = TRUE,
  tip_labels_states = FALSE,
  tip_labels_states_size = 2,
  tip_labels_states_offset = 0.1,
  node_labels_as = NULL,
  node_labels_size = 2,
  node_labels_offset = 0.1,
  pie_colors = "default",
  node_pie_size = 1,
  shoulder_pie_size = node_pie_size,
  tip_pies = TRUE,
  tip_pie_size = 0.5,
  node_pie_nudge_x = 0,
  node_pie_nudge_y = 0,
  tip_pie_nudge_x = node_pie_nudge_x,
  tip_pie_nudge_y = node_pie_nudge_y,
  shoulder_pie_nudge_x = node_pie_nudge_x,
  shoulder_pie_nudge_y = node_pie_nudge_y,
  state_transparency = 0.75,
  timeline = FALSE,
  geo = timeline,
  geo_units = list("epochs", "periods"),
  time_bars = timeline,
  ...
)

Arguments

Value

A ggplot object

Examples

# Standard ancestral state reconstruction example

# process file and assign state labels
file <- system.file("extdata",
                    "comp_method_disc/ase_freeK.tree",
                    package="RevGadgets")
example <- processAncStates(file,
                            state_labels = c("1" = "Awesome",
                                             "2" = "Beautiful",
                                             "3" = "Cool!"))
# plot (this may take a while)
plotAncStatesPie(t = example)

# DEC Biogeographic range evolution example (with timeline)

# process file
file <- system.file("extdata", "dec/simple.ase.tre", package="RevGadgets")

# labels that correspond to each region/ possible combination of regions
labs <- c("1" = "K", "2" = "O", "3" = "M", "4" = "H", "5" = "KO",
          "6" = "KM", "7" = "OM", "8" = "KH", "9" = "OH", "10" = "MH",
          "11" = "KOM", "12" = "KOH", "13" = "KMH", "14" = "OMH",
          "15" = "KOMH")
dec_example  <- processAncStates(file , state_labels = labs)
# Use the state_labels in returned tidytree object to define color palette
# These state_labels may be a subset of the labels you provided
# (not all possible regions may be sampled in the dataset)
colors <- colorRampPalette(colFun(12))(length(dec_example@state_labels))
names(colors) <- dec_example@state_labels

# plot
plotAncStatesPie(t = dec_example, pie_colors = colors, tip_labels_size = 3,
        cladogenetic = TRUE, tip_labels_offset = 0.25, timeline = TRUE,
        geo = FALSE) +
        ggplot2::theme(legend.position = c(0.1, 0.75))

Description

Plots the probability distribution of the number of lineages through time inferred with the Occurrence Birth Death Process #'

Usage

plotDiversityOBDP(
  Kt_mean,
  xlab = "Time",
  ylab = "Number of lineages",
  xticks_n_breaks = 5,
  col_Hidden = "dodgerblue3",
  col_LTT = "gray25",
  col_Total = "forestgreen",
  col_Hidden_interval = "dodgerblue2",
  col_Total_interval = "darkolivegreen4",
  palette_Hidden = c("transparent", "dodgerblue2", "dodgerblue3", "dodgerblue4", "black"),
  palette_Total = c("transparent", "green4", "forestgreen", "black"),
  line_size = 0.7,
  interval_line_size = 0.5,
  show_Hidden = TRUE,
  show_LTT = TRUE,
  show_Total = TRUE,
  show_intervals = TRUE,
  show_densities = TRUE,
  show_expectations = TRUE,
  use_interpolate = TRUE
)

Arguments

Value

A ggplot object

Examples

## Not run: 
# first run readOBDP()
start_time_trace_file <- 
     system.file("extdata", "obdp/start_time_trace.p", package="RevGadgets")
popSize_distribution_matrices_file <- 
     system.file("extdata", "obdp/Kt_trace.p", package="RevGadgets")
trees_trace_file <- 
     system.file("extdata", "obdp/mcmc_OBDP_trees.p", package="RevGadgets")
    
Kt_mean <- readOBDP( start_time_trace_file=start_time_trace_file, 
                     popSize_distribution_matrices_file=popSize_distribution_matrices_file, 
                     trees_trace_file=trees_trace_file )

# then get the customized ggplot object with plotDiversityOBDP()
p <- plotDiversityOBDP( Kt_mean,
                        xlab="Time (My)",
                        ylab="Number of lineages",
                        xticks_n_breaks=21,
                        col_Hidden="dodgerblue3",
                        col_LTT="gray25",
                        col_Total="forestgreen",
                        col_Hidden_interval="dodgerblue2",
                        col_Total_interval="darkolivegreen4",
                        palette_Hidden=c("transparent", "dodgerblue2", "dodgerblue3", 
                                         "dodgerblue4", "black"),
                        palette_Total=c("transparent", "green4", "forestgreen", "black"),
                        line_size=0.7,
                        interval_line_size=0.5,
                        show_Hidden=TRUE,
                        show_LTT=TRUE,
                        show_Total=TRUE,
                        show_intervals=TRUE,
                        show_densities=TRUE,
                        show_expectations=TRUE,
                        use_interpolate=TRUE )

# basic plot
p

# option: add a stratigraphic scale
library(deeptime)
library(ggplot2)
q <- gggeo_scale(p, dat="periods", height=unit(1.3, "line"), abbrv=F, size=4.5, neg=T)
r <- gggeo_scale(q, dat="epochs", height=unit(1.1, "line"), abbrv=F, size=3.5, neg=T, 
                    skip=c("Paleocene", "Pliocene", "Pleistocene", "Holocene"))
s <- gggeo_scale(r, dat="stages", height=unit(1, "line"), abbrv=T, size=2.5, neg=T)
s

## End(Not run)

Description

Plots the output of an episodic diversification rate analysis

Usage

plotDivRates(rates, facet = TRUE)

Arguments

Details

Plots the output of episodic diversification rate analyses. Takes as input the output of processDivRates() and plotting parameters. For now, only variable names (under "item") that contain the word "rate" are included in the plot.

The return object can be manipulated. For example, you can change the axis labels, the color palette, whether the axes are to be linked, or the overall plotting style/theme, just as with any ggplot object.

Value

A ggplot object

Examples

# download the example datasets to working directory

url_ex_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_times.log"
dest_path_ex_times <- "primates_EBD_extinction_times.log"
download.file(url_ex_times, dest_path_ex_times)

url_ex_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_rates.log"
dest_path_ex_rates <- "primates_EBD_extinction_rates.log"
download.file(url_ex_rates, dest_path_ex_rates)

url_sp_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_times.log"
dest_path_sp_times <- "primates_EBD_speciation_times.log"
download.file(url_sp_times, dest_path_sp_times)

url_sp_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_rates.log"
dest_path_sp_rates <- "primates_EBD_speciation_rates.log"
download.file(url_sp_rates, dest_path_sp_rates)

# to run on your own data, change this to the path to your data file
speciation_time_file <- dest_path_sp_times
speciation_rate_file <- dest_path_sp_rates
extinction_time_file <- dest_path_ex_times
extinction_rate_file <- dest_path_ex_rates

rates <- processDivRates(speciation_time_log = speciation_time_file,
                         speciation_rate_log = speciation_rate_file,
                         extinction_time_log = extinction_time_file,
                         extinction_rate_log = extinction_rate_file,
                         burnin = 0.25)

# then plot results:
p <- plotDivRates(rates = rates);p

# change the x-axis
p <- p + ggplot2::xlab("Thousands of years ago");p

# change the colors
p <- p + ggplot2::scale_fill_manual(values = c("red",
                                               "green",
                                               "yellow",
                                               "purple")) +
  ggplot2::scale_color_manual(values = c("red",
                                         "green",
                                         "yellow",
                                         "purple"));p

# let's say we don't want to plot relative-extinction rate,
# and use the same y-axis for all three rates
rates <- rates[!grepl("relative-extinction", rates$item),]
p2 <- plotDivRates(rates)
p2 <- p2 + ggplot2::facet_wrap(ggplot2::vars(item), scale = "fixed");p2

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_sp_times, dest_path_ex_times,
            dest_path_sp_rates, dest_path_ex_rates)

Description

Plots a single FBD tree, such as an MCC or MAP tree.

Usage

plotFBDTree(
  tree,
  timeline = FALSE,
  geo = timeline,
  geo_units = list("epochs", "periods"),
  time_bars = timeline,
  node_age_bars = TRUE,
  tip_age_bars = TRUE,
  age_bars_color = "blue",
  age_bars_colored_by = NULL,
  age_bars_width = 1.5,
  node_labels = NULL,
  node_labels_color = "black",
  node_labels_size = 3,
  node_labels_offset = 0,
  tip_labels = TRUE,
  tip_labels_italics = FALSE,
  tip_labels_formatted = FALSE,
  tip_labels_remove_underscore = TRUE,
  tip_labels_color = "black",
  tip_labels_size = 3,
  tip_labels_offset = 0,
  node_pp = FALSE,
  node_pp_shape = 16,
  node_pp_color = "black",
  node_pp_size = "variable",
  branch_color = "black",
  color_branch_by = NULL,
  line_width = 1,
  label_sampled_ancs = FALSE,
  ...
)

Arguments

Details

Plots a single tree, such as an MCC or MAP tree, with special features for plotting the output of fossilized birth-death analyses.

Value

returns a single plot object.

Examples

file <- system.file("extdata", "fbd/bears.mcc.tre", package="RevGadgets")
tree <- readTrees(paths = file)
plotFBDTree(tree = tree, timeline = TRUE, tip_labels_italics = FALSE,
            tip_labels_remove_underscore = TRUE,
            node_age_bars = TRUE, age_bars_colored_by = "posterior",
            age_bars_color = rev(colFun(2))) +
  ggplot2::theme(legend.position=c(.25, .85))

Description

plotHiSSE

Usage

plotHiSSE(rates)

Arguments

Value

a ggplot object

Examples

# download the example dataset to working directory

url <- "https://revbayes.github.io/tutorials/intro/data/primates_HiSSE_2.log"
dest_path <- "primates_HiSSE_2.log"
download.file(url, dest_path)

# to run on your own data, change this to the path to your data file
hisse_file <- dest_path

pdata <- processSSE(hisse_file)
p <- plotHiSSE(pdata);p

# change colors:
p + ggplot2::scale_fill_manual(values = c("red","green"))

# change x-axis label
p + ggplot2::xlab("Rate (events/Ma)")

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path)

Description

Plots the support (as 2ln Bayes factors) for mass extinctions.

Usage

plotMassExtinctions(
  mass_extinction_trace,
  mass_extinction_times,
  mass_extinction_name,
  prior_prob,
  return_2lnBF = TRUE
)

Arguments

Details

Works only for analyses with a fixed grid where mass extinctions may occur.

The return object can be manipulated. For example, you can change the axis labels, the color palette, whether the axes are to be linked, or the overall plotting style/theme, just as with any ggplot object.

Value

A ggplot object

References

Kass and Raftery (1995) Bayes Factors. JASA, 90 (430), 773-795.

Examples

# download the example dataset to working directory
url <-
  "https://revbayes.github.io/tutorials/intro/data/crocs_mass_extinction_probabilities.log"
dest_path <- "crocs_mass_extinction_probabilities.log"
download.file(url, dest_path)

# to run on your own data, change this to the path to your data file
mass_extinction_probability_file <- dest_path

mass_extinction_probabilities <-
              readTrace(mass_extinction_probability_file,burnin = 0.25)

# prior probability of mass extinction at any time
prior_n_expected <- 0.1
n_intervals <- 100
prior_prob <- prior_n_expected/(n_intervals-1)

# times when mass extinctions were allowed
tree_age <- 243.5
interval_times <- tree_age * seq(1/n_intervals,(n_intervals-1) /
                                                   n_intervals,1/n_intervals)

# then plot results:
p <- plotMassExtinctions(mass_extinction_trace=mass_extinction_probabilities,
                         mass_extinction_times=interval_times,
                         mass_extinction_name="mass_extinction_probabilities"
                         ,prior_prob);p

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path)

Description

plotMuSSE

Usage

plotMuSSE(rates)

Arguments

Value

a ggplot object

Examples

# download the example dataset to working directory

url <-
  "https://revbayes.github.io/tutorials/intro/data/primates_BiSSE_activity_period.log"
dest_path <- "primates_BiSSE_activity_period.log"
download.file(url, dest_path)

# to run on your own data, change this to the path to your data file
bisse_file <- dest_path

pdata <- processSSE(bisse_file)
p <- plotMuSSE(pdata);p

# change colors:
p + ggplot2::scale_fill_manual(values = c("red","green"))

# change x-axis label
p + ggplot2::xlab("Rate (events/Ma)")

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path)

Description

Plots the output of a coalescent demographic analysis.

Usage

plotPopSizes(
  df,
  plot_CIs = TRUE,
  add = FALSE,
  existing_plot = NULL,
  col = "#00883a"
)

Arguments

Details

Plots the output of coalescent demographic analyses. Takes as input the output of processPopSizes() and plotting parameters.

The return object can be manipulated. For example, you can change the axis labels, the color palette, whether the axes are to be linked, or the overall plotting style/theme, just as with any ggplot object.

Value

a ggplot object

Examples

df <- dplyr::tibble("time" = c(0.0, 1.0, 2.0, 3.0, 4.0),
                    "value" = c(1.0, 1.5, 2.0, 1.5, 1.5),
                    "upper" = c(3.5, 7.0, 6.5, 5.0, 5.0),
                    "lower" = c(0.5, 0.1, 0.5, 0.5, 0.8))

plotPopSizes(df)

Description

Plots the posterior predictive statistics data

Usage

plotPostPredStats(
  data,
  prob = c(0.9, 0.95),
  col = NULL,
  side = "both",
  type = "strict",
  PPES = FALSE,
  ...
)

Arguments

Details

Produces one ggplot object per metric. Intended to plot the results of the RevBayes tutorial: Assessing Phylogenetic Reliability Using RevBayes and P3 Model adequacy testing using posterior prediction (Data Version).

Each plot shows the rejection region for the provided quantiles, as well as a p-value for the observed statistic. If side="left" (or "right"), then the p-value is the fraction of simulated statistics that are less than ( or greater than) or equal to the observed statistic. If side="both", then the p-value is calculated by first fitting a KDE to the samples, then computing the fraction of simulated statistics with density lower than the density of he observed statistic; in this sense, the "both" option computes the size of HPD defined by the observed statistic.

Value

A list of ggplot objects, where each plot contains a density distribution of the predicted values and a dashed line of the empirical value. The blue shaded region of the density plot corresponds to the 5% two-sided quantile and the orange corresponds to the 2% two-sided quantile.

Examples

# download the example datasets to working directory

url_emp <-
   "https://revbayes.github.io/tutorials/intro/data/empirical_data_pps_example.csv"
dest_path_emp <- "empirical_data_pps_example.csv"
download.file(url_emp, dest_path_emp)

url_sim <-
   "https://revbayes.github.io/tutorials/intro/data/simulated_data_pps_example.csv"
dest_path_sim <- "simulated_data_pps_example.csv"
download.file(url_sim, dest_path_sim)

# to run on your own data, change this to the path to your data file
file_sim <- dest_path_sim
file_emp <- dest_path_emp

t <- processPostPredStats(path_sim = file_sim,
                         path_emp = file_emp)
plots <- plotPostPredStats(data = t)
plots[[1]]

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_sim, dest_path_emp)

Description

Plots the posterior distributions of variables from trace file.

Usage

plotTrace(trace, color = "default", vars = NULL, match = NULL)

Arguments

Details

Plots the posterior distributions of continuous variables from one or multiple traces (as in, from multiple runs). Shaded regions under the curve represent the 95% credible interval. If multiple traces are provided, plotTrace() will plot each run independently as well as plot the combined output. Note that for variables ith very different distributions, overlaying the plots may result in illegible figures. In these cases, we recommend plotting each parameter separately.

Value

plotTrace() returns a list of the length of provided trace object, plus one combined trace. Each element of the list contains a ggplot object with plots of the provided parameters. These plots may be modified in typical ggplot fashion.

Examples

# example with quantitative parameters

# download the example dataset to working directory
url_gtr <-
   "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.log"
dest_path_gtr <- "primates_cytb_GTR.log"
download.file(url_gtr, dest_path_gtr)

# to run on your own data, change this to the path to your data file
file <- dest_path_gtr

one_trace <- readTrace(paths = file)
plots <- plotTrace(trace = one_trace,
                    vars = c("pi[1]","pi[2]","pi[3]","pi[4]"))
plots[[1]]

# add custom colors
plots <- plotTrace(trace = one_trace,
                   vars = c("pi[3]","pi[4]","pi[1]","pi[2]"),
                   color = c("pi[1]" = "green",
                             "pi[2]"= "red",
                             "pi[3]"= "blue",
                             "pi[4]"= "orange"))
plots[[1]]

# make the same plot, using match
plots <- plotTrace(trace = one_trace, match = "pi")
plots[[1]]

#' # remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_gtr)

# plot some qualitative variables

# download the example dataset to working directory
url_rj <- "https://revbayes.github.io/tutorials/intro/data/freeK_RJ.log"
dest_path_rj <- "freeK_RJ.log"
download.file(url_rj, dest_path_rj)

file <- dest_path_rj
trace <- readTrace(path = file)

plots <- plotTrace(trace = trace,
                   vars = c("prob_rate_12", "prob_rate_13",
                            "prob_rate_31", "prob_rate_32"))
plots[[1]]

# with custom colors
plots <- plotTrace(trace = trace,
                   vars = c("prob_rate_12", "prob_rate_13",
                            "prob_rate_31", "prob_rate_32"),
                   color = c("prob_rate_12" = "green",
                             "prob_rate_13" = "red",
                             "prob_rate_31"= "blue",
                             "prob_rate_32" = "orange"))
plots[[1]]

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_rj)

Description

Plots a single tree, such as an MCC or MAP tree.

Usage

plotTree(
  tree,
  timeline = FALSE,
  geo_units = list("epochs", "periods"),
  geo = timeline,
  time_bars = timeline,
  node_age_bars = FALSE,
  age_bars_color = "blue",
  age_bars_colored_by = NULL,
  age_bars_width = 1.5,
  node_labels = NULL,
  node_labels_color = "black",
  node_labels_size = 3,
  node_labels_offset = 0,
  tip_labels = TRUE,
  tip_labels_italics = FALSE,
  tip_labels_formatted = FALSE,
  tip_labels_remove_underscore = TRUE,
  tip_labels_color = "black",
  tip_labels_size = 3,
  tip_labels_offset = 0,
  node_pp = FALSE,
  node_pp_shape = 16,
  node_pp_color = "black",
  node_pp_size = "variable",
  branch_color = "black",
  color_branch_by = NULL,
  line_width = 1,
  tree_layout = "rectangular",
  ...
)

Arguments

Details

Plots a single tree, such as an MCC or MAP tree, with optionally labeled posterior probabilities at nodes, a timescale plotted on the x - axis, and 95% CI for node ages.

Value

returns a single plot object.

Examples

# Example of standard tree plot

file <- system.file("extdata",
                    "sub_models/primates_cytb_GTR_MAP.tre",
                    package="RevGadgets")
tree <- readTrees(paths = file)
# Reroot tree before plotting
tree_rooted <- rerootPhylo(tree = tree, outgroup = "Galeopterus_variegatus")
# Plot
p <- plotTree(tree = tree_rooted, node_labels = "posterior");p

# Plot unladderized tree
p <- plotTree(tree = tree_rooted,
              node_labels = "posterior",
              ladderize = FALSE);p

# We can add a scale bar:
p + ggtree::geom_treescale(x = -0.35, y = -1)

# Example of coloring branches by rate
file <- system.file("extdata",
                    "relaxed_ou/relaxed_OU_MAP.tre",
                    package="RevGadgets")
tree <- readTrees(paths = file)
p <- plotTree(tree = tree,
              node_age_bars = FALSE,
              node_pp = FALSE,
              tip_labels_remove_underscore = TRUE,
              tip_labels_italics = FALSE,
              color_branch_by = "branch_thetas",
              line_width = 1.7) +
       ggplot2::theme(legend.position=c(.1, .9));p

Description

Plots a tree, such as an MCC or MAP tree

Usage

plotTreeFull(
  tree,
  timeline,
  geo,
  geo_units,
  time_bars,
  node_age_bars,
  tip_age_bars,
  age_bars_color,
  age_bars_colored_by,
  age_bars_width,
  node_labels,
  node_labels_color,
  node_labels_size,
  node_labels_offset,
  tip_labels,
  tip_labels_italics,
  tip_labels_formatted,
  tip_labels_remove_underscore,
  tip_labels_color,
  tip_labels_size,
  tip_labels_offset,
  label_sampled_ancs,
  node_pp,
  node_pp_shape,
  node_pp_color,
  node_pp_size,
  branch_color,
  color_branch_by,
  line_width,
  tree_layout,
  ...
)

Arguments

Details

Plots a single tree, such as an MCC or MAP tree, with the full set of functionality to be called by plotTree() and plotFBDTree()

Value

returns a single plot object.

See Also

called by plotTree and plotFBDTree

Description

Turn posterior samples collected by MCMC into a parametric prior distribution.

Usage

posteriorSamplesToParametricPrior(
  samples,
  distribution,
  variance_inflation_factor = 2
)

Arguments

Details

The distributions are fit by the method of moments. The function allows inflating the prior variance relative to the posterior being supplied.

Value

Numeric vector of parameters with names (to avoid rate/scale and var/sd confusion).

Examples

# download the example datasets to working directory

url_ex_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_times.log"
dest_path_ex_times <- "primates_EBD_extinction_times.log"
download.file(url_ex_times, dest_path_ex_times)

url_ex_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_rates.log"
dest_path_ex_rates <- "primates_EBD_extinction_rates.log"
download.file(url_ex_rates, dest_path_ex_rates)

url_sp_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_times.log"
dest_path_sp_times <- "primates_EBD_speciation_times.log"
download.file(url_sp_times, dest_path_sp_times)

url_sp_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_rates.log"
dest_path_sp_rates <- "primates_EBD_speciation_rates.log"
download.file(url_sp_rates, dest_path_sp_rates)

# to run on your own data, change this to the path to your data file
speciation_time_file <- dest_path_sp_times
speciation_rate_file <- dest_path_sp_rates
extinction_time_file <- dest_path_ex_times
extinction_rate_file <- dest_path_ex_rates

primates <- processDivRates(speciation_time_log = speciation_time_file,
                            speciation_rate_log = speciation_rate_file,
                            extinction_time_log = extinction_time_file,
                            extinction_rate_log = extinction_rate_file,
                            burnin = 0.25)

speciation_rates <-
       dplyr::pull(primates[which(primates$item == "speciation rate"),],
                  "value")
speciation_1_gamma_prior <-
          posteriorSamplesToParametricPrior(speciation_rates,"gamma")

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_sp_times, dest_path_ex_times,
            dest_path_sp_rates, dest_path_ex_rates)

Description

Process data for ancestral states plotting

Usage

processAncStates(
  path,
  state_labels = NULL,
  labels_as_numbers = FALSE,
  missing_to_NA = TRUE
)

Arguments

Value

A treedata object

Examples

# standard ancestral state estimation example
file <- system.file("extdata",
                    "comp_method_disc/ase_freeK.tree",
                    package="RevGadgets")
example <- processAncStates(file,
                            state_labels = c("1" = "Awesome",
                                             "2" = "Beautiful",
                                             "3" = "Cool!"))

#chromosome evolution example
file <- system.file("extdata",
                    "chromo/ChromEvol_simple_final.tree",
                    package="RevGadgets")
chromo_example <- processAncStates(file, labels_as_numbers = TRUE)

Description

processBranchData

Usage

processBranchData(
  tree,
  dat,
  burnin = 0.25,
  parnames = c("avg_lambda", "avg_mu", "num_shifts"),
  summary = "median",
  net_div = FALSE
)

Arguments

Value

a treedata file with attached branch-specific data

Examples

# download the example dataset to working directory
url_rates <-
  "https://revbayes.github.io/tutorials/intro/data/primates_BDS_rates.log"
dest_path_rates <- "primates_BDS_rates.log"
download.file(url_rates, dest_path_rates)

url_tree <-
  "https://revbayes.github.io/tutorials/divrate/data/primates_tree.nex"
dest_path_tree <- "primates_tree.nex"
download.file(url_tree, dest_path_tree)

# to run on your own data, change this to the path to your data file
treefile <- dest_path_tree
logfile <- dest_path_rates

branch_data <- readTrace(logfile)
tree <- readTrees(paths = treefile)

annotated_tree <- processBranchData(tree, branch_data, summary = "median")

# you can plot this output
p <- plotTree(tree = annotated_tree,
              node_age_bars = FALSE,
              node_pp = FALSE,
              tip_labels = FALSE,
              color_branch_by = "avg_lambda",
              line_width = 0.8) +
     ggplot2::theme(legend.position=c(.1, .9));p
# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_tree, dest_path_rates)

Description

Processing the output of a episodic diversification rate analysis with mass-extinction events.

Usage

processDivRates(
  speciation_time_log = "",
  speciation_rate_log = "",
  extinction_time_log = "",
  extinction_rate_log = "",
  fossilization_time_log = "",
  fossilization_rate_log = "",
  burnin = 0.25,
  probs = c(0.025, 0.975),
  summary = "median"
)

Arguments

Details

For processing the output of an episodic diversification rate analysis. processDivRates() assumes that the epochs are fixed rather than inferred. Additionally, it assumes that times correspond to rates such that the first rate parameter (i.e. speciation[1]) corresponds to the present. Conversely, the first time parameter (i.e. interval_times[1]) corresponds to the first time interval after the present, moving backwards in time. processDivRates() relies on readTrace and produces a list object that can be read by plotDivRates() to visualize the results. For now, only one log file per parameter type is accepted (i.e. log files from multiple runs must be combined before reading into the function).

Value

List object with processed rate and time parameters.

Examples

# download the example datasets to working directory

url_ex_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_times.log"
dest_path_ex_times <- "primates_EBD_extinction_times.log"
download.file(url_ex_times, dest_path_ex_times)

url_ex_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_extinction_rates.log"
dest_path_ex_rates <- "primates_EBD_extinction_rates.log"
download.file(url_ex_rates, dest_path_ex_rates)

url_sp_times <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_times.log"
dest_path_sp_times <- "primates_EBD_speciation_times.log"
download.file(url_sp_times, dest_path_sp_times)

url_sp_rates <-
   "https://revbayes.github.io/tutorials/intro/data/primates_EBD_speciation_rates.log"
dest_path_sp_rates <- "primates_EBD_speciation_rates.log"
download.file(url_sp_rates, dest_path_sp_rates)

# to run on your own data, change this to the path to your data file
speciation_time_file <- dest_path_sp_times
speciation_rate_file <- dest_path_sp_rates
extinction_time_file <- dest_path_ex_times
extinction_rate_file <- dest_path_ex_rates

rates <- processDivRates(speciation_time_log = speciation_time_file,
                         speciation_rate_log = speciation_rate_file,
                         extinction_time_log = extinction_time_file,
                         extinction_rate_log = extinction_rate_file,
                         burnin = 0.25)

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_sp_times, dest_path_ex_times,
            dest_path_sp_rates, dest_path_ex_rates)

Description

Processing the output of a coalescent demographic analysis.

Usage

processPopSizes(
  population_size_log = "",
  interval_change_points_log = "",
  model = "constant",
  burnin = 0.25,
  probs = c(0.025, 0.975),
  summary = "median",
  num_grid_points = 100,
  spacing = "exponential",
  max_age = NULL,
  min_age = NULL,
  distribution = FALSE
)

Arguments

Details

For processing the output of a coalescent demographic analysis. processPopSizes() assumes that the the first size parameter (i.e. population_size[1]) corresponds to the present. processPopSizes() partly relies on readTrace and produces a list object that can be read by plotPopSizes() to visualize the results. For now, only one log file per parameter type is accepted (i.e. log files from multiple runs must be combined before reading into the function).

Value

List object with processed rate and, if applicable, time parameters (if distribution = FALSE). Matrix object with distributions of population size (if distribution = TRUE). If applicable, one row for each point on the grid, with the times of the grid points as row names.

Description

Reads in and processes posterior-predictive statistics

Usage

processPostPredStats(path_sim, path_emp)

Arguments

Value

A list of data frames

Examples

# download the example datasets to working directory

url_emp <-
   "https://revbayes.github.io/tutorials/intro/data/empirical_data_pps_example.csv"
dest_path_emp <- "empirical_data_pps_example.csv"
download.file(url_emp, dest_path_emp)

url_sim <-
   "https://revbayes.github.io/tutorials/intro/data/simulated_data_pps_example.csv"
dest_path_sim <- "simulated_data_pps_example.csv"
download.file(url_sim, dest_path_sim)

# to run on your own data, change this to the path to your data file
file_sim <- dest_path_sim
file_emp <- dest_path_emp

t <- processPostPredStats(path_sim = file_sim,
                         path_emp = file_emp)

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_sim, dest_path_emp)

Description

Title

Usage

processSSE(
  path,
  speciation = "speciation",
  extinction = "extinction",
  speciation_hidden = "speciation_hidden",
  rates = c(speciation, extinction, "net-diversification"),
  ...
)

Arguments

Value

a data frame

Examples

# download the example dataset to working directory

url <-
  "https://revbayes.github.io/tutorials/intro/data/primates_BiSSE_activity_period.log"
dest_path <- "primates_BiSSE_activity_period.log"
download.file(url, dest_path)

# to run on your own data, change this to the path to your data file
bisse_file <- dest_path

pdata <- processSSE(bisse_file)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path)

Description

Reads and formats the outputs of an analysis with the Occurrence Birth Death Process (MCMC parameter inference + diversity estimation)

Usage

readOBDP(
  start_time_trace_file,
  popSize_distribution_matrices_file,
  trees_trace_file
)

Arguments

Value

A data.frame

Examples

## Not run: 
# first run readOBDP()
start_time_trace_file <- 
     system.file("extdata", "obdp/start_time_trace.p", package="RevGadgets")
popSize_distribution_matrices_file <- 
     system.file("extdata", "obdp/Kt_trace.p", package="RevGadgets")
trees_trace_file <- 
     system.file("extdata", "obdp/mcmc_OBDP_trees.p", package="RevGadgets")
    
Kt_mean <- readOBDP( start_time_trace_file=start_time_trace_file, 
                     popSize_distribution_matrices_file=popSize_distribution_matrices_file, 
                     trees_trace_file=trees_trace_file )

# then get the customized ggplot object with plotDiversityOBDP()
p <- plotDiversityOBDP( Kt_mean,
                        xlab="Time (My)",
                        ylab="Number of lineages",
                        xticks_n_breaks=21,
                        col_Hidden="dodgerblue3",
                        col_LTT="gray25",
                        col_Total="forestgreen",
                        col_Hidden_interval="dodgerblue2",
                        col_Total_interval="darkolivegreen4",
                        palette_Hidden=c("transparent", "dodgerblue2", "dodgerblue3", 
                                         "dodgerblue4", "black"),
                        palette_Total=c("transparent", "green4", "forestgreen", "black"),
                        line_size=0.7,
                        interval_line_size=0.5,
                        show_Hidden=TRUE,
                        show_LTT=TRUE,
                        show_Total=TRUE,
                        show_intervals=TRUE,
                        show_densities=TRUE,
                        show_expectations=TRUE,
                        use_interpolate=TRUE )

# basic plot
p

# option: add a stratigraphic scale
library(deeptime)
library(ggplot2)
q <- gggeo_scale(p, dat="periods", height=unit(1.3, "line"), abbrv=F, size=4.5, neg=T)
r <- gggeo_scale(q, dat="epochs", height=unit(1.1, "line"), abbrv=F, size=3.5, neg=T, 
                    skip=c("Paleocene", "Pliocene", "Pleistocene", "Holocene"))
s <- gggeo_scale(r, dat="stages", height=unit(1, "line"), abbrv=T, size=2.5, neg=T)
s

## End(Not run)

Description

Reads in MCMC log files

Usage

readTrace(
  paths,
  format = "simple",
  delim = "\t",
  burnin = 0.1,
  check.names = FALSE,
  ...
)

Arguments

Details

Reads in one or multiple MCMC log files from the same analysis and discards a user-specified burn-in, compatible with multiple monitor types. If the trace contains vectors of vectors and the user does not specify format = "complex", readTrace() will read in those columns as factors rather than as numeric vectors.

Value

List of dataframes (of length 1 if only 1 log file provided).

Examples

# read and process a single trace file


# download the example dataset to working directory
url_gtr <-
   "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.log"
dest_path_gtr <- "primates_cytb_GTR.log"
download.file(url_gtr, dest_path_gtr)

# to run on your own data, change this to the path to your data file
file_single <- dest_path_gtr

one_trace <- readTrace(paths = file_single)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_gtr)

# read and process multiple trace files, such as from multiple runs of
# the same analysis

# download the example dataset to working directory
url_1 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_1.log"
dest_path_1 <- "primates_cytb_GTR_run_1.log"
download.file(url_1, dest_path_1)

url_2 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_2.log"
dest_path_2 <- "primates_cytb_GTR_run_2.log"
download.file(url_2, dest_path_2)

# to run on your own data, change this to the path to your data file
file_1 <- dest_path_1
file_2 <- dest_path_2

# read in the multiple trace files
multi_trace <- readTrace(path = c(file_1, file_2), burnin = 0.0)

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_1, dest_path_2)

Description

Reads in a tree file containing one or multiple trees

Usage

readTrees(paths, tree_name = "psi", burnin = 0, n_cores = 1L, verbose = TRUE)

Arguments

Details

Reads in a tree file in either nexus or newick format, and containing a single tree or multiple trees (as in the results of a Bayesian analysis). For reading in annotated tree files of continuous character evolution, the parameter must be considered a node parameter rather than branch parameter. Set isNodeParameter = TRUE in the extended newick monitor (mnExtNewick)

Value

A list (across runs) of lists (across samples) of treedata objects.

Examples

# read in a single nexus file

# download the example dataset to working directory
url_nex <-
 "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_MAP.tre"
dest_path_nex <- "primates_cytb_GTR_MAP.tre"
download.file(url_nex, dest_path_nex)

# to run on your own data, change this to the path to your data file
file <- dest_path_nex
tree_single_old <- readTrees(paths = file)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_nex)

# read in a single newick string

# download the example dataset to working directory
url_new <-
 "https://revbayes.github.io/tutorials/intro/data/primates.tre"
dest_path_new <- "primates.tre"
download.file(url_new, dest_path_new)

# to run on your own data, change this to the path to your data file
file_new <- dest_path_new
tree_new <- readTrees(paths = file_new)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_new)


# read in a tree trace (may take a few seconds)

# download the example dataset to working directory
url_multi <-
 "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.trees"
dest_path_multi <- "primates_cytb_GTR.trees"
download.file(url_multi, dest_path_multi)

# to run on your own data, change this to the path to your data file
file_multi <- dest_path_multi
tree_multi <- readTrees(paths = file_multi)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_multi)

Description

Removes burnin from MCMC trace

Usage

removeBurnin(trace, burnin)

Arguments

Details

Removes burnin from an MCMC trace, such as the output of readTrace(). If multiple traces are provided, this function will remove the burnin from each.

Value

List of dataframes (of length 1 if only 1 log file provided).

Examples

# download the example dataset to working directory
url_gtr <-
    "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.log"
dest_path_gtr <- "primates_cytb_GTR.log"
download.file(url_gtr, dest_path_gtr)

# to run on your own data, change this to the path to your data file
file_single <- dest_path_gtr

one_trace <- readTrace(paths = file_single)
one_trace_burnin <- removeBurnin(trace = one_trace, burnin = 0.1)

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_gtr)

Description

Reroots a phylogeny given an outgroup taxon or clade

Usage

rerootPhylo(tree, outgroup)

Arguments

Details

Modifies a tree object by rerooting using a specified outgroup taxon or clade. Places the root at the midpoint of the branch subtending the outgroup. If the input contains multiple trees, all trees will be rerooted.

Value

returns a list of list of treedata objects, with the trees rooted.

See Also

phytools: reroot.

Examples

file <- system.file("extdata",
                    "sub_models/primates_cytb_GTR_MAP.tre",
                    package="RevGadgets")
tree <- readTrees(paths = file)
# root with one taxon
tree_rooted <- rerootPhylo(tree = tree, outgroup = "Galeopterus_variegatus")
# root with clade, specified by two taxa
tree_rooted <- rerootPhylo(tree = tree,
                           outgroup = c("Varecia_variegata_variegata",
                                        "Propithecus_coquereli"))

Description

This package provides functions to process and plot the output of RevBayes analyses.

Description

This function finds the global scale parameter value that produces the desired prior mean number of "effective" rate shifts. Given a specified magnitude for an effective shift, shift_size, an effective shift occurs when two adjacent values are more than shift_size-fold apart from each other. That is, an effective shift is the event that rate[i+1]/rate[i] > shift_size or rate[i+1]/rate[i] < 1/shift_size.

Usage

setMRFGlobalScaleHyperpriorNShifts(
  n_episodes,
  model,
  prior_n_shifts = log(2),
  shift_size = 2
)

Arguments

Details

Finding these values for a HSMRF model can take several seconds for large values of n_episodes because of the required numerical integration.

Value

The hyperprior.

References

Magee et al. (2019) Locally adaptive Bayesian birth-death model successfully detects slow and rapid rate shifts. doi: https://doi.org/10.1101/853960

Examples

# Get global scale for a HSMRF model with 100 episodes.
gs <- setMRFGlobalScaleHyperpriorNShifts(100, "HSMRF")

# Plot a draw from this HSMRF distribution

trajectory <- simulateMRF(n_episodes = 100,
                          model = "HSMRF",
                          global_scale_hyperprior = gs)

plot(1:100,
     rev(trajectory),
     type = "l",
     xlab = "time",
     ylab = "speciation rate")

Description

This function simulates a draw from a HSMRF or GMRF distribution given a user-specified global scale parameter. The MRF can be taken to be on the log-scale (such as for a birth rate) or the real-scale. The first value must be specified

Usage

simulateMRF(
  n_episodes,
  model,
  global_scale_hyperprior,
  initial_value = NULL,
  exponentiate = TRUE
)

Arguments

Value

A vector drawn from the specified MRF model on the specified (log- or real-) scale.

References

Magee et al. (2020) Locally adaptive Bayesian birth-death model successfully detects slow and rapid rate shifts. PLoS Computational Biology, 16 (10): e1007999.

Faulkner, James R., and Vladimir N. Minin. Locally adaptive smoothing with Markov random fields and shrinkage priors. Bayesian analysis, 13 (1), 225.

Examples

# Simulate a 100-episode HSMRF model for a speciation-rate through time
trajectory <- simulateMRF(n_episodes = 100,
                          model = "HSMRF",
                          global_scale_hyperprior = 0.0021)
plot(1:100,
     rev(trajectory),
     type = "l",
     xlab = "time",
     ylab = "speciation rate")

Description

Summarizes trace file(s) that have been read into memory

Usage

summarizeTrace(trace, vars)

Arguments

Details

Summarizes a trace file for continuous or discrete characters by computing the mean and 95% credible interval for quantitative character and the 95% credible set for discrete characters.

Value

summarizeTrace() returns a list of the length of provided variables. For quantitative variables, it returns the mean and 95 For discrete variables, it returns the 95 associated probabilities.

Examples

# continuous character only example, one run

# download the example dataset to working directory
url_gtr <-
    "https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR.log"
dest_path_gtr <- "primates_cytb_GTR.log"
download.file(url_gtr, dest_path_gtr)

# to run on your own data, change this to the path to your data file
file_single <- dest_path_gtr

one_trace <- readTrace(paths = file_single)
trace_sum <- summarizeTrace(trace = one_trace,
                            vars = c("pi[1]","pi[2]","pi[3]","pi[4]"))
trace_sum[["pi[1]"]]

# remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_gtr)

# continuous character example, multiple runs

#' # download the example dataset to working directory
url_1 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_1.log"
dest_path_1 <- "primates_cytb_GTR_run_1.log"
download.file(url_1, dest_path_1)

url_2 <-
"https://revbayes.github.io/tutorials/intro/data/primates_cytb_GTR_run_2.log"
dest_path_2 <- "primates_cytb_GTR_run_2.log"
download.file(url_2, dest_path_2)

# to run on your own data, change this to the path to your data file
file_1 <- dest_path_1
file_2 <- dest_path_2

# read in the multiple trace files
multi_trace <- readTrace(path = c(file_1, file_2), burnin = 0.0)

trace_sum_multi <- summarizeTrace(trace = multi_trace,
                                  vars = c("pi[1]","pi[2]","pi[3]","pi[4]"))
trace_sum_multi[["pi[1]"]]

# remove files
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_1, dest_path_2)

# discrete character example

# download the example dataset to working directory
url_rj <- "https://revbayes.github.io/tutorials/intro/data/freeK_RJ.log"
dest_path_rj <- "freeK_RJ.log"
download.file(url_rj, dest_path_rj)

file <- dest_path_rj
trace <- readTrace(path = file)

trace_sum_discrete <- summarizeTrace(trace = trace,
                                     vars = c("prob_rate_12",
                                              "prob_rate_13",
                                              "prob_rate_31",
                                              "prob_rate_32"))
trace_sum_discrete[["prob_rate_12"]]

#' # remove file
# WARNING: only run for example dataset!
# otherwise you might delete your data!
file.remove(dest_path_rj)