makePrior: Generate prior distributions for the multivariate Brownian...

In Caetanods/ratematrix: Bayesian Estimation of the Evolutionary Rate Matrix

Generate prior distributions for the multivariate Brownian motion model

Description

Generates prior densities for the MCMC sampler.

Usage

makePrior(
  r,
  p,
  den.mu = "unif",
  par.mu,
  den.sd = "unif",
  par.sd,
  unif.corr = TRUE,
  Sigma = NULL,
  nu = NULL
)

Arguments

r	number of traits in the model.
p	number of evolutionary rate matrix regimes fitted to the phylogeny.
den.mu	one of "unif" (uniform prior, default) or "norm" (normal prior).
par.mu	the parameters for the prior density on the vector of phylogenetic means. Matrix with 2 columns and number of rows equal to the number of traits (r). When the density ('den.mu') is set to "unif" then par.mu[,1] is the minimum values and par.mu[,2] is the maximum values for each trait. When the density is set to "norm" then par.mu[,1] is the mean values and par.mu[,2] is the standard deviation values for the set of normal densities around the vector of phylogenetic means.
den.sd	one of "unif" (uniform prior, default) or "lnorm" (log-normal prior).
par.sd	the parameters for the density of standard deviations. Matrix with 2 columns and number of rows equal to the number of evolutionary rate matrix regimes fitted to the phylogenetic tree (p). When "den.sd" is set to "unif", then 'par.sd[,1]' (the minimum) need to be a vector of positive values and 'par.sd[,2]' is the vector of maximum values. When "den.sd" is set to "lnorm" then 'par.sd[,1]' is the vector of log(means) for the density and 'par.sd[,2]' is the vector of log(standard deviations) for the distributions. If there is only one regime fitted to the tree, then 'par.sd' is a vector with length 2 (e.g., c(min, max)).
unif.corr	whether the correlation structure of the prior distribution on the Sigma matrix is flat. This sets an uniformative prior (as uninformative as possible) to the evolutionary correlations among the traits.
Sigma	the scale matrix to be used as a parameter for the inverse-Wishart distribution responsible for the generation of the correlation matrices. Need to be a list of matrices with number of elements equal to the number of evolutionary rate regimes fitted to the phylogeny, in the case of a single regime, this needs to be a matrix (not a list). The scale matrix is somewhat analogous to the mean of a normal distribution. Thus, if 'Sigma' show strong positive correlation among traits, then the distribution generated by the prior will vary around positive correlations. This parameter will be ignored if 'unif.corr' is set to TRUE.
nu	the degrees of freedom parameter for the inverse-Wishart distribution. Need to be a numeric vector with length equal to the number of evolutionary rate matrix regimes fitted to tree. Larger values of 'nu' will make the prior distribution closer around 'Sigma' and small values will increase the variance. This parameter is analogous to the variance parameter of a normal distribution, however, it has an inverted relationship.

Details

This function is integrated within the 'ratematrixMCMC' function that runs the MCMC chain. However, this implementation allows for more control over the prior distribution for the analysis. The prior functions produced here can be easily passed to the 'ratematrixMCMC' function. See examples and more information in the 'ratematrixMCMC' function.

One can use the output of this function in order to sample from the prior using the 'samplePrior' function. A sample from the prior can be set as the starting point of the MCMC sampler.

Independent priors are defined for the phylogenetic mean, the vector of standard deviations and the structure of correlation, allowing for a wide range of configurations. Priors for the phylogenetic mean and the standard deviations can be uniform or normal (lognormal in the case of the standard deviations). The prior on the matrix of correlations is distributed as an inverse-Wishart and can be set to a marginally uniform prior or to be centered around a given variance-covariace matrix.

The prior for the model has two elements, one is the vector of phylogenetic means (or the root values) and the other is the evolutionary rate matrices (the vcv matrices for the rate of the multivariate BM model). The vector of root values can be distributed as any continuous distribution. In this implementation the two options are the uniform and the normal distribution. On the other hand, the prior distribution for the rate matrices need to be more elaborated. Here we divide the variance-covariance matrix into two elements, a correlation matrix and a vector of standard deviations. Standard deviations can be modelled as any continuous distrbuted of positive values. Here we use a uniform or a log-normal distribution. The correlation matrix need to be derived from a distribution of covariance matrices known as the inverse-Wishart. The inverse-Wishart is controlled by two parameters; the scale matrix (Sigma) and the degrees of freedom (nu). Any variance-covariance matrix can be used as the scale matrix. To set a marginally uniform prior for the correlation structure of the evolutionary rate matrices sampled for the model one need to set 'Sigma' as an identity matrix and 'nu' as the dimension of the matrix +1. This is performed automatically by the function when the option 'unif.corr' is set to TRUE.

Value

List of density functions to compute the log prior probability of parameter values.

Author(s)

Daniel S. Caetano and Luke J. Harmon

Examples

data( centrarchidae )
## Set the limits of the uniform prior on the root based on the observed traits
data.range <- t( apply( centrarchidae$data, 2, range ) )
## Set a reasonable value for the uniform prior distribution for the standard deviation.
## Here the minimum rate for the traits is 0 and the maximum is 10 ( using 'sqrt(10)' to
##      transform to standard deviation).
## Note that we need to use a matrix with dimension dependent on the number of traits.
par.sd <- cbind(c(0,0), sqrt( c(10,10) ))
prior <- makePrior(r = 2, p = 2, den.mu = "unif", par.mu = data.range, den.sd = "unif"
                   , par.sd = par.sd)
## Running a very short chain, it will not converge:
handle <- ratematrixMCMC(data=centrarchidae$data, phy=centrarchidae$phy.map, prior=prior
                         , gen=5000, dir=tempdir())
posterior <- readMCMC(handle, burn = 0.2, thin = 1)
plotRatematrix( posterior )

Caetanods/ratematrix documentation built on Jan. 4, 2023, 3:12 p.m.