ABCSMC {SimBIID}

Runs ABC-SMC algorithm

Description

Runs the Approximate Bayesian Computation Sequential Monte Carlo (ABC-SMC) algorithm of Toni et al. (2009) for fitting infectious disease models to time series count data.

Usage

ABCSMC(x, ...)

## S3 method for class 'ABCSMC'
ABCSMC(
  x,
  tols = NULL,
  ptols = NULL,
  mintols = NULL,
  ngen = 1,
  parallel = FALSE,
  mc.cores = NA,
  ...
)

## Default S3 method:
ABCSMC(
  x,
  priors,
  func,
  u,
  tols = NULL,
  ptols = NULL,
  mintols = NULL,
  ngen = 1,
  npart = 100,
  parallel = FALSE,
  mc.cores = NA,
  ...
)

Arguments

x

An ABCSMC object or a named vector with entries containing the observed summary statistics to match to. Names must match to 'tols'.

...

Further arguments to pass to func. (Not used if extending runs.)

tols

A vector or matrix of tolerances, with the number of rows defining the number of generations required, and columns defining the summary statistics to match to. If a vector, then the length determines the summary statistics. The columns/entries must match to those in 'x'.

ptols

The proportion of simulated outcomes at each generation to use to derive adaptive tolerances.

mintols

A vector of minimum tolerance levels.

ngen

The number of generations of ABC-SMC to run.

parallel

A logical determining whether to use parallel processing or not.

mc.cores

Number of cores to use if using parallel processing.

priors

A data.frame containing columns parnames, dist, p1 and p2, with number of rows equal to the number of parameters. The column parname simply gives names to each parameter for plotting and summarising. Each entry in the dist column must contain one of c("unif", "norm", "gamma"), and the corresponding p1 and p2 entries relate to the hyperparameters (lower and upper bounds in the uniform case; mean and standard deviation in the normal case; and shape and rate in the gamma case).

func

Function that runs the simulator and checks whether the simulation matches the data. The first four arguments must be pars, data, tols and u. If the simulations do not match the data then the function must return an NA, else it must returns a vector of simulated summary measures. In this latter case the output from the function must be a vector with length equal to ncol(data) and with entries in the same order as the columns of data.

u

A named vector of initial states.

npart

An integer specifying the number of particles.

Details

Samples initial particles from the specified prior distributions and then runs a series of generations of ABC-SMC. The generations can either be specified with a set of fixed tolerances, or by setting the tolerances at each new generation as a quantile of the tolerances of the accepted particles at the previous generation. Uses bisection method as detailed in McKinley et al. (2018). Passing an ABCSMC object into the ABCSMC() function acts as a continuation method, allowing further generations to be run.

Value

An ABCSMC object, essentially a list containing:

• pars: a list of matrix objects containing the accepted particles. Each element of the list corresponds to a generation of ABC-SMC, with each matrix being of dimension npart x npars;

• output: a list of matrix objects containing the simulated summary statistics. Each element of the list corresponds to a generation of ABC-SMC, with each matrix being of dimension npart x ncol(data);

• weights: a list of vector objects containing the particle weights. Each element of the list corresponds to a generation of ABC-SMC, with each vector being of length npart;

• ESS: a list of effective sample sizes. Each element of the list corresponds to a generation of ABC-SMC, with each vector being of length npart;

• accrate: a vector of length nrow(tols) containing the acceptance rates for each generation of ABC;

• tols: a copy of the tols input;

• ptols: a copy of the ptols input;

• mintols: a copy of the mintols input;

• priors: a copy of the priors input;

• data: a copy of the data input;

• func: a copy of the func input;

• u a copy of the u input;

• addargs: a copy of the ... inputs.

References

Toni T, Welch D, Strelkowa N, Ipsen A and Stumpf MP (2009) <doi:10.1098/rsif.2008.0172>

McKinley TJ, Cook AR and Deardon R (2009) <doi:10.2202/1557-4679.1171>

McKinley TJ, Vernon I, Andrianakis I, McCreesh N, Oakley JE, Nsubuga RN, Goldstein M and White RG (2018) <doi:10.1214/17-STS618>

See Also

print.ABCSMC, plot.ABCSMC, summary.ABCSMC

Examples

## set up SIR simulationmodel
transitions <- c(
    "S -> beta * S * I -> I", 
    "I -> gamma * I -> R"
)
compartments <- c("S", "I", "R")
pars <- c("beta", "gamma")
model <- mparseRcpp(
    transitions = transitions, 
    compartments = compartments,
    pars = pars
)
model <- compileRcpp(model)

## generate function to run simulators
## and return summary statistics
simSIR <- function(pars, data, tols, u, model) {

    ## run model
    sims <- model(pars, 0, data[2] + tols[2], u)
    
    ## this returns a vector of the form:
    ## completed (1/0), t, S, I, R (here)
    if(sims[1] == 0) {
        ## if simulation rejected
        return(NA)
    } else {
        ## extract finaltime and finalsize
        finaltime <- sims[2]
        finalsize <- sims[5]
    }
    
    ## return vector if match, else return NA
    if(all(abs(c(finalsize, finaltime) - data) <= tols)){
        return(c(finalsize, finaltime))
    } else {
        return(NA)
    }
}

## set priors
priors <- data.frame(
    parnames = c("beta", "gamma"), 
    dist = rep("gamma", 2), 
    stringsAsFactors = FALSE
)
priors$p1 <- c(10, 10)
priors$p2 <- c(10^4, 10^2)

## define the targeted summary statistics
data <- c(
    finalsize = 30, 
    finaltime = 76
)

## set initial states (1 initial infection 
## in population of 120)
iniStates <- c(S = 119, I = 1, R = 0)

## set initial tolerances
tols <- c(
    finalsize = 50,
    finaltime = 50
)

## run 2 generations of ABC-SMC
## setting tolerance to be 50th
## percentile of the accepted 
## tolerances at each generation
post <- ABCSMC(
    x = data, 
    priors = priors, 
    func = simSIR, 
    u = iniStates, 
    tols = tols, 
    ptol = 0.2, 
    ngen = 2, 
    npart = 50,
    model = model
)
post

## run one further generation
post <- ABCSMC(post, ptols = 0.5, ngen = 1)
post
summary(post)

## plot posteriors
plot(post)

## plot outputs
plot(post, "output")