run.evorates: Run sampler for an Evolving Rates model

In bstaggmartin/backwards-BM-simulator: Evolving Rate Models of Continuous Trait Evolution

Run sampler for an Evolving Rates model

Description

This function runs a Stan-based Hamiltonian Monte Carlo (HMC) sampler given input data and priors, but only stores, rather than processes, the output.

Usage

run.evorates(
  input.evorates.obj,
  return.as.obj = TRUE,
  out.file = NULL,
  check.overwrite = TRUE,
  constrain.Rsig2 = FALSE,
  trend = FALSE,
  lik.power = 1,
  ...
)

Arguments

input.evorates.obj	The output of input.evorates()
return.as.obj	TRUEor FALSE: should results be passed back to the R environment as a evorates_fit object? Defaults to TRUE. If FALSE, results are saved to files with automatically generated names instead (see below).
out.file	A directory to save results to. If unspecified, an automatic directory is generated based on the current date and time and R's working directory. The function will generate csv files for each chain of the HMC sampler using Stan's built-in functionality (note that these will thus be on the transformed scale), as well as a separate RDS file giving additional information about the data. Defaults to NULL, which means no results are saved to file, though this will be changed automatically if return.as.obj = FALSE.
check.overwrite	TRUE or FALSE: should files in the directory specified by out.file be checked to prevent accidentally overwriting existing files? This part of the code is not thoroughly tested and might take a long time for folders with many files, so some users may wish to just switch it off. Defaults to TRUE.
constrain.Rsig2	TRUE or FALSE: should the rate varinace (R_sig2) parameter be constrained to 0, resulting in a simple Brownian Motion or early/late burst model? Defaults to FALSE. See Details for a definition of model parameters.
trend	TRUE or FALSE: should a trend in rates over time (R_mu) be estimated, resulting in an early/late burst or trended evorates model? Defaults to FALSE.
lik.power	A single number between 0 and 1 specifying what power to raise the likelihood function to. This is useful for sampling from "power posteriors" that shift the posterior to look more like the prior (indeed, you can set this to 0 to sample from the prior itself). Useful for model diagnostics and calculating things like Bayes Factors. Technically, you can set lik.power above 1 (a technique called "data cloning") but this is not beneficial for evorates models and we don't recommend it.
...	Other optional arguments: Prior arguments (see details for further information on what parameters mean): In general, prior standard deviations and, in some cases, means, can be tweaked by passing arguments named "<parameter_name>_sd" and "<parameter_name>_mean", respectively. For example: R_mu_mean = 1 or R_sig2_sd = 1. Note that all prior arguments are necessarily on the sampling scale when passed to run.evorates() Prior on rate variance (R_sig2): Follows a half-Cauchy distribution (basically a half-normal distribution with extremely fat tails) with adjustable standard deviation. By default, standard deviation is set to 5 divided by the maximum height of tree (i.e., 5 on the sampling scale). Prior on trend (R_mu): Follows a normal distribution with adjustable mean and standard deviation. By default, mean is set to 0 (also 0 on sampling scale), and standard deviation is set to 10 divided by the maximum height of tree (i.e., 10 on the sampling scale). Prior on rate at the root (R_0): Follows a normal distribution with adjustable mean and standard deviation. By default, mean is set to variance of trait.data divided by the maximum height of tree on the natural log scale (i.e., 0 on sampling scale), and standard deviation is set to 10 (also 10 on sampling scale). Prior on tip error variance (Y_sig2): Follows a truncated t distribution (basically a normal distribution with potentially fatter tails) with adjustable mean, standard deviation, and degrees of freedom. The distribution is truncated at 0 to the left such that Y_sig2 estimates are always positive (note that setting the mean to 0 results in a half-t distribution). By default, mean is set to 0, standard deviation is set to 2 times the variance of trait.data (i.e., 2 on the sampling scale), and degrees of freedom is set to 1. The degrees of freedom can be tweaked with the argument "Ysig2_df" and may be any positive number, including Inf, with lower numbers leading to a less informative priors with fatter tails. Notably, setting "Ysig2_df" to Inf makes this prior a truncated normal and much more informative! Additional arguments to pass to rstan::sampling(), most commonly: chains to specify the number of HMC chains (defaults to 4) iter to specify the number of iterations in HMC chains (defaults to 2000) warmup to specify the number of warmup iterations (defaults to floor(iter/2)) thin to specify which iterations to keep in results (defaults to 1 or no thinning) cores to specify the number of computer cores to use (defaults to getOption("mc.cores", 1L)) refresh to control when progress is reported (defaults to max(iter/10, 1), and can be suppressed by setting to 0 or less) There are other things users might want to mess with, like seed, init, and control. See ?rstan::sampling and ?rstan::stan for further details.

Value

A list if return.as.obj = TRUE, containing the unprocessed model fit in stanfit format and updated input.evorates.obj. Otherwise, the directory results were saved to (see out.file for details).

#' @details Parameter definitions:

• Rate variance (R_sig2, also denoted \sigma^2_{\sigma^2}): Determines how much random variation accumulates in rates over time. Specifically, independent lineages evolving for a length of time t would exhibit log-normally distributed rates with standard deviation sqrt(t * R_sig2). Another way to think about this is that the 95% credible interval for fold-changes in rates over 1 unit of time is given by exp(c(-1,1) * 1.96 * sqrt(R_sig2)), assuming R_mu = 0.

• Trend (R_mu, also denoted \mu_{\sigma^2}): Determines whether median rates tend to decrease (if negative) or increase (if positive) over time. Specifically, the median fold-change in rate for a given lineage is exp(t * R_mu) over a length of time t. This is distinct from changes in average rates, which depends on both R_mu and R_sig2.

• Note on combining rate variance and trend parameters: The 95% credible interval for fold-changes in rates over 1 unit of time, given a trend, is simply exp(R_mu) * exp(c(-1,1) * 1.96 * sqrt(R_sig2)) or equivalently exp(R_mu + c(-1,1) * 1.96 * sqrt(R_sig2)). This distribution of rate change is right-skewed due to the exp() function. Because of this, median changes in rates over time will always be lower than average changes. The "average change" parameter, denoted R_del or \delta_{\sigma^2}, is given by R_mu + R_sig2 / 2. Accordingly, the average fold-change in rate for a given lineage is exp(t * R_del) over a length of time t.

• Tip error variance (Y_sig2, also denoted \sigma^2_y): Determines the among of "error" around observed trait values for tips without fixed standard errors. Specifically, raw observations for a given tip are sampled from a normal distribution centered at that tip's "true trait value" with standard deviation sqrt(Y_sig2).

• Rate at the root (R_0, also denoted \sigma^2_0): The natural log of the rate at the root of the entire phylogeny.

• Branchwise rates (R_i, also denoted ln \bar{\sigma^2_i}, where i is the index of an edge in tree): The natural log of the average rate along branches of the phylogeny. Note that rates are always shifting over time under this model, and these quantities are thus averages. The true rate value at any particular time point along a branch is a related, but separate, quantity. These will be NA for branches of length 0.

• Background rate (bg_rate): The natural log of the average trait evolution rate for the entire phylogeny, given by the average of exp(R_i), with each entry weighted by its respective branch length. NAs are ignored in this calculation since they correspond to branches of length 0.

• Branchwise rate deviations (Rdev_i, also denoted ln \bar{\sigma^2_{dev,i}}): Determines whether a branches exhibit relatively "fast" (if positive) or "slow" rates (if negative). Generally, these are the differences between the branchwise rates and geometric background rate on the natural log scale, though this depends on remove.trend). Here, the geometric background rate is defined as the weighted average of R_i, with weights corresponding to branch lengths. This helps prevent some issues with comparing highly right-skewed distributions. If remove.trend = TRUE, then branchwise rates are first "detrended" prior to calculating background rates and deviations (R_i - (-log(abs(R_mu)) - log(l_i) + log(abs(exp(R_mu * t1_i) - exp(R_mu * t0_i)))), where l is a vector of branch lengths and t0/t1 are vectors of start and end times of each branch). This basically make branchwise rate deviations determine whether branches exhibit slow/fast rates given the overall trend in rates through time. Otherwise, branchwise rate deviations will simmply indicate slow/fast branches occur at the root/tips of a tree in the presence of a strong trend.

See Also

fit.evorates is a convenient wrapper for input.evorates, run.evorates, and output.evorates

Examples

#get whale/dolphin tree/trait data
data("cet_fit")
tree <- cet_fit$call$tree
trait.data <- cet_fit$call$trait.data

#prepare data for evorates model
input <- input.evorates(tree = tree, trait.data = trait.data)
#run the evorates model sampler (takes a couple minutes)
run <- run.evorates(input.evorates.obj = input, out.file = "test_fit", chains = 1)
#process the output from the object
fit <- output.evorates(run.evorates.obj = run)
#process the output from the files
fit <- output.evorates(run.evorates.obj = "test_fit")

#see fit.evorates() documentation for further examples

bstaggmartin/backwards-BM-simulator documentation built on June 3, 2024, 5:51 p.m.