fit.evorates: Fit an Evolving Rates model

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

Fit an Evolving Rates model

Description

This function processes tree and trait data, runs a Stan-based Hamiltonian Monte Carlo (HMC) sampler to fit these data to an evorates model, and returns the output of this sampler in a (relatively) user-friendly format.

Usage

fit.evorates(
  tree,
  trait.data,
  trait.se = NULL,
  constrain.Rsig2 = FALSE,
  trend = FALSE,
  lik.power = 1,
  sampling.scale = FALSE,
  return.as.obj = TRUE,
  out.file = NULL,
  check.overwrite = TRUE,
  include.warmup = FALSE,
  report.quantiles = c(0.025, 0.5, 0.975),
  report.means = TRUE,
  report.devs = TRUE,
  remove.trend = TRUE,
  ...
)

Arguments

tree	An object of class "phylo"
trait.data	Three options: A named vector of trait values. A rownamed matrix of trait values. For the moment, must be 1 column since multivariate models are not yet supported. A data.frame with 2 columns: 1 with numeric data which is interpreted as trait values and 1 with string/factor data interpreted as names. If more than 1 of either kind of column is found, returns error. Will also use rownames if no string/factor column is found, but this limits data to consist of only 0-1 observations per tip. Both multiple observations for a single tip and missing observations are allowed. In all cases, the associated names must generally match the tip labels found in tree (tree$tip.label) exactly. The only exception are numbers corresponding to internal node indices (these can be viewed by running: plot(tree); nodelabels()), which are used to assign trait values to internal nodes. As of now, the function actually completely ignores tree$node.label–I should probably fix that and plan to in the future. Be aware that this feature slightly alters tree because it works by grafting extra tips of length 0 to the appropriate nodes.
trait.se	A vector, matrix, or data.frame of trait value standard errors which must be unambiguously labeled (see trait.data). Alternatively, a single, unlabeled number that will be applied to all tips in tree (e.g., you could set it to 0 to specify all trait values are known without error). If NULL (the default), standard error is estimated for all tips in the tree while fitting the model. There are several things to note here: Unlike trait.data, there can only be 1 trait standard error per tip. You can use NA to specify tips that you want to estimate standard error for and Inf to specify tips that have missing trait values. Generally, tips with multiple observations and no observations are automatically assigned standard errors of NA and Inf, respectively. Any tips with unspecified standard error default to NA. Any conflicting/impossible standard error specifications are corrected and return warnings.
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.
sampling.scale	TRUE or FALSE: should provided prior parameters (see ...) be interpreted on the raw scale of the data or on the transformed scale? All data passed to the Stan-based HMC sampler are transformed such that tree's total height is 1 and the variance of the trait data is 1. Defaults to FALSE, such that prior parameters are interpreted on the untransformed scale.
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.
include.warmup	TRUE or FALSE: should warmup be included in posterior samples? Warmup is always included for parameters used to tune the HMC chain, but warmup may be included or excluded for actual estimated parameters. Defaults to FALSE.
report.quantiles	A vector of posterior distribution quantiles to return (should be between 0 and 1). Set to NULL to not return any quantiles. Defaults to 2.5%, 50%, and 97.5% quantiles.
report.means	TRUE or FALSE: should posterior distribution means be returned? Defaults to TRUE.
report.devs	TRUE or FALSE: should the difference between branchwise rates and the overall "background rate" on the natural log scale ("rate deviations") be returned? The background rate is simply the mean branchwise rate, weighted by each branch's respective length. If TRUE, also calculates the posterior probability rate deviations are greater than 0. These additional parameters help give a sense of which branches in tree exhibit anomalous trait evolution rates. Defaults to TRUE, but is automatically switched to FALSE when fitting a Brownian Motion model (constrain.Rsig2 = FALSE & trend = FALSE).
remove.trend	TRUE or FALSE: should the rate deviation calculations remove (i.e., "account for") the effect of trends in rates over time? This may be helpful since strong trends in rates can mask otherwise anomalously high or low trait evolution rates in certain parts of the tree. Defaults to TRUE, but has no effect if no trend was fitted or if report.devs is FALSE. If TRUE, note that report.devs will be switched to FALSE if fitting early/late burst model (constrain.Rsig2 = FALSE & trend = TRUE).
...	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. 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.

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.

Value

An object of class "evorates_fit" if return.as.obj = TRUE. Otherwise, the directory results were saved to (see out.file for details). An evorates_fit object is a list of at least 5 components:

• call, which contains information on the final tree, trait values, trait standard errors, and prior parameters passed to Stan's HMC sampling algorithm (on untransformed scale for better interpretability, see sampling.scale).

• sampler.control, which contains various information on the HMC run, including the number of chains, iterations, warmup, thinning rate, etc.

• sampler.params, an array of class "param_block", containing parameters/diagnostics that were used to tune the behavior of the HMC while Stan ran it, as well as the (log) prior (prior__), likelihood (lik__), and posterior probability (post__) of each iteration in the HMC. See Stan manual for more information on what the parameters mean. The likelihood is not raised to lik.power here, but the log posterior probability is calculated while accounting for lik.power. This always includes warmup iterations, though these can be discarded using exclude.warmup() or combine.chains().

• param.diags, a param_block array, containing diagnostics for each parameter estimated during the fit, including the initial value of the HMC chain (init), the bulk effective sample size (bulk_ess), the tail effective sample size (tail_ess), and the Rhat (Rhat). See Rhat for more details on what these diagnostics mean. Generally, you want effective sample sizes to be greater than 100 times the number of chains and Rhat to be less than 1.01. The functions check.ess() and check.mix() will check these thresholds for you automatically.

• chains, another param_block array containing sampled parameter values for each parameter estimated during the fit. See details for further information on what each parameter means.

• The object optionally contains more param_block arrays of posterior distribution quantiles (quantiles) and means (means), and posterior probabilities (post.prob, see report.devs and remove.trend).

All param_block arrays' dimensions go in the order of iterations/diagnostics/quantiles, then parameters, then chains.

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

#fit data to evorates model (takes a couple minutes)
fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1)

#specifying parameter constraints
#add trend parameter
trend.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1,
                          trend = TRUE)
#only trend parameter (early burst model)
EB.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1,
                       trend = TRUE, constrain.Rsig2 = TRUE)
#no Rsig2 or trend (Brownian Motion model)
BM.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1,
                       constrain.Rsig2 = TRUE)

#specifying trait standard error
#by default, estimates standard error, but you can set all tips to have no standard error
noSE.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1,
                         trait.se = 0)
#or even set it for specific tips
SE.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1
                       trait.se = setNames(c(1, 2), c("Orcinus_orca", "Balaenoptera_musculus")))
                       
#specifying priors
#make Rsig2 prior tighter (normally it's 5)
Risg2.fit <- fit.evorates(tree  = tree, trait.data = trait.data, chains = 1,
                          R_sig2_prior_sd = 1, sampling.scale = TRUE)
#maybe you REALLY believe rates should decrease by ~1% every million years?
Rmu.fit <- fit.evorates(tree = tree, trait.data = trait.data, chains = 1,
                        trend = TRUE,
                        R_mu_prior_mean = log(1 - 0.01), R_mu_prior_sd = log(1 + 0.005))

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