PMCMC: Runs particle MCMC algorithm
In SimBIID: Simulation-Based Inference Methods for Infectious Disease Models

Description Usage Arguments Details Value References See Also Examples

Runs particle Markov chain Monte Carlo (PMCMC) algorithm using a bootstrap particle filter for fitting infectious disease models to time series count data.

PMCMC(x, ...)

## S3 method for class 'PMCMC'
PMCMC(
  x,
  niter = 1000,
  nprintsum = 100,
  adapt = TRUE,
  adaptmixprop = 0.05,
  nupdate = 100,
  ...
)

## Default S3 method:
PMCMC(
  x,
  priors,
  func,
  u,
  npart = 100,
  iniPars = NA,
  fixpars = FALSE,
  niter = 1000,
  nprintsum = 100,
  adapt = TRUE,
  propVar = NA,
  adaptmixprop = 0.05,
  nupdate = 100,
  ...
)

`x`	A `PMCMC` object, or a `data.frame` containing time series count data, with the first column called `t`, followed by columns of time-series counts. The time-series counts columns must be in the order of the 'counts' object in the 'func' function (see below).
`...`	Not used here.
`niter`	An integer specifying the number of iterations to run the MCMC.
`nprintsum`	Prints summary of MCMC to screen every `nprintsum` iterations. Also defines how often adaptive scaling of proposal variances occur.
`adapt`	Logical determining whether to use adaptive proposal or not.
`adaptmixprop`	Mixing proportion for adaptive proposal.
`nupdate`	Controls when to start adaptive update.
`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`	A `SimBIID_model` object or an `XPtr` to simulation function. If the latter, then this function must take the following arguments in order: `NumericVector pars`: a vector of parameters; `double tstart`: the start time; `double tstop`: the end time; `IntegerVector u`: a vector of states at time `tstart`; `IntegerVector counts`: a vector of observed counts at `tstop`.
`u`	A named vector of initial states.
`npart`	An integer specifying the number of particles for the bootstrap particle filter.
`iniPars`	A named vector of initial values for the parameters of the model. If left unspecified, then these are sampled from the prior distribution(s).
`fixpars`	A logical determining whether to fix the input parameters (useful for determining the variance of the marginal likelihood estimates).
`propVar`	A numeric (npars x npars) matrix with log (or logistic) covariances to use as (initial) proposal matrix. If left unspecified then defaults to `diag(nrow(priors)) * (0.1 ^ 2) / nrow(priors)`.

Function runs a particle MCMC algorithm using a bootstrap particle filter for a given model. If running with fixpars = TRUE then this runs niter simulations using fixed parameter values. This can be used to optimise the number of particles after a training run. Also has print(), summary(), plot(), predict() and window() methods.

If the code throws an error, then it returns a missing value (NA). If fixpars = TRUE it returns a list of length 2 containing:

output: a matrix with two columns. The first contains the simulated log-likelihood, and the second is a binary indicator relating to whether the simulation was 'skipped' or not (1 = skipped, 0 = not skipped);
pars: a vector of parameters used for the simulations.

If fixpars = FALSE, the routine returns a PMCMC object, essentially a list containing:

pars: an mcmc object containing posterior samples for the parameters;
u: a copy of the u input;
accrate: the cumulative acceptance rate;
npart: the chosen number of particles;
time: the time taken to run the routine (in seconds);
propVar: the proposal covariance for the parameter updates;
data: a copy of the x input;
priors: a copy of the priors input;
func: a copy of the func input.

Andrieu C, Doucet A and Holenstein R (2010) <doi:10.1111/j.1467-9868.2009.00736.x>

print.PMCMC, plot.PMCMC, predict.PMCMC, summary.PMCMC window.PMCMC

## set up data to pass to PMCMC
flu_dat <- data.frame(
    t = 1:14,
    Robs = c(3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29, 14, 4)
)

## set up observation process
obs <- data.frame(
    dataNames = "Robs",
    dist = "pois",
    p1 = "R + 1e-5",
    p2 = NA,
    stringsAsFactors = FALSE
)

## set up model (no need to specify tspan
## argument as it is set in PMCMC())
transitions <- c(
    "S -> beta * S * I / (S + I + R + R1) -> I", 
    "I -> gamma * I -> R",
    "R -> gamma1 * R -> R1"
)
compartments <- c("S", "I", "R", "R1")
pars <- c("beta", "gamma", "gamma1")
model <- mparseRcpp(
    transitions = transitions, 
    compartments = compartments,
    pars = pars,
    obsProcess = obs
)

## set priors
priors <- data.frame(
    parnames = c("beta", "gamma", "gamma1"), 
    dist = rep("unif", 3), 
    stringsAsFactors = FALSE)
priors$p1 <- c(0, 0, 0)
priors$p2 <- c(5, 5, 5)

## define initial states
iniStates <- c(S = 762, I = 1, R = 0, R1 = 0)

set.seed(50)

## run PMCMC algorithm
post <- PMCMC(
    x = flu_dat, 
    priors = priors,
    func = model, 
    u = iniStates,
    npart = 25,
    niter = 5000, 
    nprintsum = 1000
)

## plot MCMC traces
plot(post, "trace")

## continue for some more iterations
post <- PMCMC(post, niter = 5000, nprintsum = 1000)

## plot traces and posteriors
plot(post, "trace")
plot(post)

## remove burn-in
post <- window(post, start = 5000)

## summarise posteriors
summary(post)