A comment on a recent post gave me the motivation to try out the new EasyABC package for R, developed by Franck Jabot, Thierry Faure, Nicolas Dumoulin and maintained by Nicolas Dumoulin.

Approximate Bayesian Computation (ABC) is a relatively new method that allows treating any stochastic model (IBM, stochastic population model, …) in a statistical framework by generating “approximate” likelihood values through simulating from the model. We provide a gentle introduction to ABC and some alternative approaches in our recent Ecology Letters review on “statisitical inference for stochastic simulation models”.

The EasyABC package, available from CRAN, implements a number of algorithms for the three main sampling strategies used in ABC, namely Rejection Sampling, Sequential Monte Carlo (SMC) and Markov Chain Monte Carlo (MCMC). All those are also discussed in our review. The use of the package is relatively straightforward. To give a demonstration, I implemented the parameter inference of a normal distribution using the ABC-MCMC algorithm proposed by Marjoram that I coded by hand in my previous post on ABC in EasyABC. The EasyABC solution is provided below.

I think that this is an interesting package that will hopefully contribute to making ABC more widely known and used in ecology. As we say in our review, ABC has a huge potential as a solution for many typical ecological problems, but to make this more widely known is currently hindered by the fact that you have to code everything by hand, which excludes a large number of users. Thanks to the authors for their work!

library(EasyABC)

# assuming the data are 10 samples of a normal distribution
# with mean 5.3 and sd 2.7
data =  rnorm(10, mean =5.3, sd = 2.7)


# we want to use ABC to infer the parameters that were used.
# we sample from the same model and use mean and variance
# as summary statstitics for the model and the data.

# observed summary statistics
summarydata = c(mean(data), sd(data))

model <- function(par){

  # stochastic model generates a sample for given par
  samples <- rnorm(10, mean =par[1], sd = par[2])

  # returning simulated summary statistics
  return(c(mean(samples), sd(samples)))
}

# normalization of the summary statistics, 1,1 probably not the most appropriate choice but to keep it easy
tabnormalization=c(1,1)

# definition of the (flat priors)
priormatrix=cbind(c(-15,2),c(15,4))

# call EasyABC
ABC_Marjoram_original<-ABC_mcmc(method="Marjoram_original", model=model, prior_matrix=priormatrix,
    n_obs=10000, n_between_sampling=1, summary_stat_target=summarydata, dist_max=1, proposal_range=c(1,1),
tab_normalization = tabnormalization, use_seed = F)


str(ABC_Marjoram_original)
par(mfrow=c(2,1))
hist(ABC_Marjoram_original$param[5000:10000,1], main = "Posterior for slope")
hist(ABC_Marjoram_original$param[5000:10000,2], main = "Posterior for intercept")