abctools (version 0.2-2)

stage2: stage2

Description

Summary statistics selection for ABC inference using estimated posterior error.

Usage

stage2(obs, param, sumstats, obspar = NULL, chosen, dsets = 100, 
sumsubs = 1:ncol(sumstats), limit = length(sumsubs), do.only = NULL, do.err = FALSE, 
final.dens = FALSE, verbose = TRUE, ...)

Value

A list with the following components:bestthe best subset of statistics.closestthe indices of the dsets simulated datasets closest to the oberved dataset as measured by the chosen subset of summaries.errsimulation error (if obspar is supplied and do.err=TRUE).orderthe subsets considered during the algorithm (same as the input do.only.post.samplean array of dimension nacc x npar x ndatasets giving the posterior sample for each observed dataset. Not returned if final.dens=FALSE.sumsubsan index into the subsets considered during the algorithm.

Warning

This function is computationally intensive due to its cyclic ABC inference procedure.

Details

The function uses the chosen set of summaries to determine the dsets simulated datasets which are closest (in Euclidean norm) to the observed dataset. Since the model parameters generating the summary statistics are known for these simulated datasets, for each candidate subset of summary statistics, we can compute the error under ABC inference for each of these datasets. The best subset of summary statistics is that which minimizes this (average) error over all dsets datasets.

References

Blum, M. G. B, Nunes, M. A., Prangle, D. and Sisson, S. A. (2013) A comparative review of dimension reduction methods in approximate Bayesian computation. Stat. Sci. (to appear). Nunes, M. A. and Balding, D. J. (2010) On Optimal Selection of Summary Statistics for Approximate Bayesian Computation. Stat. Appl. Gen. Mol. Biol. 9, Iss. 1, Art. 34.

Examples

# load example data:

data(coal)
data(coalobs)

param<-coal[,2]
simstats<-coal[,5:8]

# use matrix below just in case to preserve dimensions.

obsstats<-matrix(coalobs[1,5:8],nrow=1)
obsparam<-matrix(coalobs[1,1])

tmp<-stage2(obsstats, param, simstats, chosen=c(1,3), dsets = 10) 
tmp$best