abc_smc_nltt {nLTT}

A function to perform Approximate Bayesian Computation within a Sequential Markov Chain (ABC-SMC), for diversification analysis of phylogenetic trees.

Description

This function performs ABC-SMC as described in Toni 2009 for given diversification model, provided a phylogenetic tree. ABC-SMC is not limited to only using the normalized LTT as statistic.

Usage

abc_smc_nltt(
    tree, statistics, simulation_function, init_epsilon_values,
    prior_generating_function, prior_density_function,
    number_of_particles = 1000, sigma = 0.05, stop_rate = 1e-05
  )

Arguments

tree

an object of class "phylo"; the tree upon which we want to fit our diversification model

statistics

A vector containing functions that take a tree as an argument and return a single scalar value (the statistic).

simulation_function

A function that implements the diversification model and returns an object of class "phylo".

init_epsilon_values

A vector containing the initial threshold values for the summary statistics from the vector statistics.

prior_generating_function

Function to generate parameters from the prior distribution of these parameters (e.g. a function returning lambda and mu in case of the birth-death model)

prior_density_function

Function to calculate the prior probability of a set of parameters.

number_of_particles

Number of particles to be used per iteration of the ABC-SMC algorithm.

sigma

Standard deviation of the perturbance distribution (perturbance distribution is a gaussian with mean 0).

stop_rate

If the acceptance rate drops below stopRate, stop the ABC-SMC algorithm and assume convergence.

Value

A matrix with n columns, where n is the number of parameters you are trying to estimate.

Author(s)

Thijs Janzen

References

Toni, T., Welch, D., Strelkowa, N., Ipsen, A., & Stumpf, M.P.H. (2009). Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface, 6(31), 187-202.

Examples

  ## Not run: 

  prior_gen <- function() {
    return( rexp(n=2, rate=0.1) )
  }

  prior_dens <- function(val) {
    return( dexp( val[1], rate = 0.1) * dexp( val[2], rate = 0.1) )
  }

  require(TESS)

  treeSim <- function(params) {
    t <- TESS.sim.age(n=1, lambda = params[1], mu = params[2], age = 10)[[1]]
    return(t)
  }

  obs <- treeSim(c(0.5,0.1))

  statWrapper <- function(tree1) {
    return( nLTTstat_exact(tree1, obs, "abs"))
  }

  stats <- c(statWrapper)

  results <-  abc.smc.nltt(
    obs, stats, treeSim, init_epsilon_values = 0.2,
    prior_generating_function = prior_gen,
    prior_density_function = prior_dens,
    number_of_particles = 1000, sigma = 0.05, stop_rate = 1e-5
  )

  
## End(Not run) # end of dontrun