tess.mcmc: tess.mcmc: Markov chain Monte Carlo simulation using a...

In hoehna/TESS: Diversification Rate Estimation and Fast Simulation of Reconstructed Phylogenetic Trees under Tree-Wide Time-Heterogeneous Birth-Death Processes Including Mass-Extinction Events

Description

tess.mcmc constructs a Markov chain Monte Carlo sampler (MCMC) by implementing a general Metropolis-Hastings algorithm. Any model can be used where the likelihood is known and thus can be passed in as an argument. The parameters have to be continuous. Proposals are taken from a normal distribution centered around the current value. The varaince of the new proposed values is initalized with 1 but can be automatically optimized when using the option adaptive = TRUE. The algorithm creates sampels from the posterior probility distribution and returns these a CODA mcmc object. More information can be obtained in the vignette about how to apply this method.

Usage

tess.mcmc(likelihoodFunction,priors,parameters,logTransforms,delta,
             iterations,burnin=round(iterations/3),thinning=1,
             adaptive=TRUE,verbose=FALSE)

Arguments

likelihoodFunction	The log-likelihood function which will be called internally by likelihoodFunction(parameters).
priors	A list of functions of the log-prior-densities of each parameter.
parameters	The initial parameter value list.
logTransforms	A vector of booleans telling if log-transform for the parameters should be used (e.g. for rates).
delta	The variance of new proposed values.
iterations	The number of iterations for the MCMC.
burnin	The number of iterations to burn before starting the MCMC.
thinning	The frequency of taking a sample of the parameters.
adaptive	Should we use adaptive MCMC?
verbose	Do you want detailed information during the run?

Value

Returns the posterior samples for the parameters.

Author(s)

Sebastian Hoehna

References

S. Hoehna: Fast simulation of reconstructed phylogenies under global, time-dependent birth-death processes. 2013, Bioinformatics, 29:1367-1374.

S. Hoehna, MR May and BR Moore: TESS: Bayesian inference of lineage diversification rates from (incompletely sampled) molecular phylogenies in R. 2015, Bioinformatics.

Examples

# load in a test data set
data(cettiidae)

# convert the phylogeny into the branching times
times <- as.numeric( branching.times(cettiidae) )

# specify a likelihood function that takes in a vector of parameters
likelihood <- function(params) {
  # We use the parameters as diversification rate and turnover rate.
  # Thus we need to transform first
  b <- params[1] + params[2]
  d <- params[2]
  
  lnl <- tess.likelihood(times,b,d,samplingProbability=1.0,log=TRUE)
  return (lnl)
}

# specify a the prior functions
prior.diversification <- function(x) { dexp(x,rate=0.1,log=TRUE) }
prior.turnover <- function(x) { dexp(x,rate=0.1,log=TRUE) }
priors <- c(prior.diversification,prior.turnover)

# Note, the number of iterations and the burnin is too small here
# and should be adapted for real analyses
samples <- tess.mcmc( likelihood,
		      priors,
		      runif(2,0,1),
		      logTransforms=c(TRUE,TRUE),
		      delta=c(0.1,0.1),
		      iterations=100,
		      burnin=20)

# now summarize and visualize the results
#plot(samples)
summary(samples)
colMeans(samples)

hoehna/TESS documentation built on Feb. 3, 2022, 5:59 a.m.