PMCMC {SimBIID}

Runs particle MCMC algorithm

Description

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

Usage

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,
  ...
)

Arguments

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).

Details

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.

Value

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.

References

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

See Also

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

Examples

## 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)