Description Usage Arguments Details Value References See Also Examples
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | 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,
...
)
|
x |
An |
... |
Further arguments to pass to |
tols |
A |
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 |
mc.cores |
Number of cores to use if using parallel processing. |
priors |
A |
func |
Function that runs the simulator and checks whether the simulation matches the data.
The first four arguments must be |
u |
A named vector of initial states. |
npart |
An integer specifying the number of particles. |
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.
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.
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>
print.ABCSMC
, plot.ABCSMC
, summary.ABCSMC
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | ## 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")
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.