abc_smc_nltt: A function to perform Approximate Bayesian Computation within...

In richelbilderbeek/nLTT: Calculate the NLTT Statistic

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

richelbilderbeek/nLTT documentation built on Aug. 23, 2023, 8 a.m.