tess.steppingStoneSampling: tess.steppingStoneSampling: Marginal likelihood estimation...

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.steppingStoneSampling uses a power posterior series and stepping-stone-sampling to estimate the marginal likelihood of a model.

Usage

tess.steppingStoneSampling(likelihoodFunction,priors,parameters,logTransforms,
                           iterations,burnin=round(iterations/3),K=50,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).
iterations	The number of iterations for the MCMC.
burnin	The number of iterations to burn before starting the MCMC.
K	The number of stepping stones.
verbose	Whether to print diagnostic outputs.

Value

Returns the posterior samples for the parameters.

Author(s)

Sebastian Hoehna

References

Xie et al., 2011: Improving marginal likelihood estimation for Bayesian phylogenetic model selection Baele et al., 2012: Improving the accuracy of demographic and molecular clock model comparison while accommodating phylogenetic uncertainty Baele et al., 2013: Accurate Model Selection of Relaxed Molecular Clocks in Bayesian Phylogenetics

Examples

data(cettiidae)
times <- as.numeric( branching.times(cettiidae) )

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)
}

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, the burnin
# and the number of stepping stones is too small here
# and should be adapted for real analyses
marginalLikelihood <- tess.steppingStoneSampling( likelihood,
						  priors,
						  runif(2,0,1),
						  c(TRUE,TRUE),
						  10,
						  10,
						  K=4)

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