pmcmc: Run a Particle MCMC Sampler within the Squire Framework
In mrc-ide/squire: SEIR transmission model of COVID-19

pmcmc

R Documentation

Run a Particle MCMC Sampler within the Squire Framework

Description

The user inputs initial parameter values for R0, Meff, and the start date The log prior likelihood of these parameters is calculated based on the user-defined prior distributions. The log likelihood of the data given the initial parameters is estimated using a particle filter, which has two functions: - Firstly, to generate a set of 'n_particles' samples of the model state space, at time points corresponding to the data, one of which is selected randomly to serve as the proposed state sequence sample at the final data time point. - Secondly, to produce an unbiased estimate of the likelihood of the data given the proposed parameters. The log posterior of the initial parameters given the data is then estimated by adding the log prior and log likelihood estimate.

The pMCMC sampler then proceeds as follows, for n_mcmc iterations: At each loop iteration the pMCMC sampler pefsorms three steps: 1. Propose new candidate samples for R0, Meff, Meff_pl, and start_date based on the current samples, using the proposal distribution (currently multivariate Gaussian with user-input covariance matrix (proposal_kernel), and reflecting boundaries defined by pars_min, pars_max) 2. Calculate the log prior of the proposed parameters, Use the particle filter to estimate log likelihood of the data given the proposed parameters, as described above, as well as proposing a model state space. Add the log prior and log likelihood estimate to estimate the log posterior of the proposed parameters given the data. 3. Metropolis-Hastings step: The joint canditate sample (consisting of the proposed parameters and state space) is then accepted with probability min(1, a), where the acceptance ratio is simply the ratio of the posterior likelihood of the proposed parameters to the posterior likelihood of the current parameters. Note that by choosing symmetric proposal distributions by including reflecting boundaries, we avoid the the need to include the proposal likelihood in the MH ratio.

If the proposed parameters and states are accepted then we update the current parameters and states to match the proposal, otherwise the previous parameters/states are retained for the next iteration.

After generating the pMCMC simulation, there are replicates specific iterations sampled based on the posterior probability. The parameters from those iterations are then used to generate new trajectories within the squire_model framework.

Usage

pmcmc(
  data,
  n_mcmc,
  log_likelihood = NULL,
  log_prior = NULL,
  n_particles = 100,
  steps_per_day = 4,
  output_proposals = FALSE,
  n_chains = 1,
  squire_model = explicit_model(),
  pars_obs = list(phi_cases = 1, k_cases = 2, phi_death = 1, k_death = 2, exp_noise =
    1e+06),
  pars_init = list(start_date = as.Date("2020-02-07"), R0 = 2.5, Meff = 2, Meff_pl = 3,
    R0_pl_shift = 0),
  pars_min = list(start_date = as.Date("2020-02-01"), R0 = 0, Meff = 1, Meff_pl = 2,
    R0_pl_shift = -2),
  pars_max = list(start_date = as.Date("2020-02-20"), R0 = 5, Meff = 3, Meff_pl = 4,
    R0_pl_shift = 5),
  pars_discrete = list(start_date = TRUE, R0 = FALSE, Meff = FALSE, Meff_pl = FALSE,
    R0_pl_shift = FALSE),
  proposal_kernel = NULL,
  scaling_factor = 1,
  reporting_fraction = 1,
  treated_deaths_only = FALSE,
  country = NULL,
  population = NULL,
  contact_matrix_set = NULL,
  baseline_contact_matrix = NULL,
  date_contact_matrix_set_change = NULL,
  R0_change = NULL,
  date_R0_change = NULL,
  hosp_bed_capacity = NULL,
  baseline_hosp_bed_capacity = NULL,
  date_hosp_bed_capacity_change = NULL,
  ICU_bed_capacity = NULL,
  baseline_ICU_bed_capacity = NULL,
  date_ICU_bed_capacity_change = NULL,
  date_vaccine_change = NULL,
  baseline_max_vaccine = NULL,
  max_vaccine = NULL,
  date_vaccine_efficacy_infection_change = NULL,
  baseline_vaccine_efficacy_infection = NULL,
  vaccine_efficacy_infection = NULL,
  date_vaccine_efficacy_disease_change = NULL,
  baseline_vaccine_efficacy_disease = NULL,
  vaccine_efficacy_disease = NULL,
  Rt_args = NULL,
  burnin = 0,
  replicates = 100,
  forecast = 0,
  required_acceptance_ratio = 0.23,
  start_adaptation = round(n_mcmc/2),
  gibbs_sampling = FALSE,
  gibbs_days = NULL,
  ...
)

Arguments

`data`	Data to fit to. This must be constructed with `particle_filter_data`
`n_mcmc`	number of mcmc mcmc iterations to perform
`log_likelihood`	function to calculate log likelihood, must take named parameter vector as input, allow passing of implicit arguments corresponding to the main function arguments. Returns a named list, with entries: - $log_likelihood, a single numeric - $sample_state, a numeric vector corresponding to the state of a single particle, chosen at random, at the final time point for which we have data. If NULL, calculated using the function calc_loglikelihood.
`log_prior`	function to calculate log prior, must take named parameter vector as input, returns a single numeric. If NULL, uses uninformative priors which do not affect the posterior
`n_particles`	Number of particles (considered for both the PMCMC fit and sampling from posterior)
`steps_per_day`	Number of steps per day
`output_proposals`	Logical indicating whether proposed parameter jumps should be output along with results
`n_chains`	number of MCMC chains to run
`squire_model`	A squire model to use
`pars_obs`	list of parameters to use in comparison with `compare`. Must be a list containing, e.g. list(phi_cases = 0.1, k_cases = 2, phi_death = 1, k_death = 2, exp_noise = 1e6)
`pars_init`	named list of initial inputs for parameters being sampled
`pars_min`	named list of lower reflecting boundaries for parameter proposals
`pars_max`	named list of upper reflecting boundaries for parameter proposals
`pars_discrete`	named list of logicals, indicating if proposed jump should be discrete
`proposal_kernel`	named matrix of proposal covariance for parameters
`scaling_factor`	numeric for starting scaling factor for covariance matrix. Default = 1
`reporting_fraction`	Reporting fraction. Numeric for what proportion of the total deaths the reported deaths represent. E.g. 0.5 results in the model calibrating to twice the deaths provided by `data$deaths`
`treated_deaths_only`	Boolean for whether likelihood is based only on deaths that occur from healthcare systems, i.e. are treated. Default = FALSE, which uses all deaths.
`country`	Character for country beign simulated. WIll be used to generate `population` and `contact_matrix_set` if unprovided. Either `country` or `population` and `contact_matrix_set` must be provided.
`population`	Population vector (for each age group). Default = NULL, which will cause population to be sourced from `country`
`contact_matrix_set`	List of contact matrices to be used from the dates provided in `date_contact_matrix_set_change`.Must be same length as `date_contact_matrix_set_change`
`baseline_contact_matrix`	The starting contact matrix prior to any changes due to interventions or otherwise. Default = NULL, which will use the contact matrix associated with the coutnry provided.
`date_contact_matrix_set_change`	Calendar dates at which the contact matrices set in `model_params` change. Defaut = NULL, i.e. no change
`R0_change`	Numeric vector for relative changes in R0. Default = NULL, i.e. no change in R0
`date_R0_change`	Calendar dates at which R0_change occurs. Defaut = NULL, i.e. no change in R0
`hosp_bed_capacity`	Number of hospital beds at each date specified in `date_hosp_bed_capacity_change`. Must be same length as `date_hosp_bed_capacity_change`.
`baseline_hosp_bed_capacity`	The starting number of hospital beds before the epidemic started. Default = NULL, which will use the hospital beds data for the country provided. If no country is provided then this is 5/1000 of the population
`date_hosp_bed_capacity_change`	Calendar dates at which hospital bed capacity changes set in `model_params` change. Defaut = NULL, i.e. no change
`ICU_bed_capacity`	Number of ICU beds at each date specified in `date_ICU_bed_capacity_change`. Must be same length as `date_ICU_bed_capacity_change`.
`baseline_ICU_bed_capacity`	The starting number of ICU beds before the epidemic started. Default = NULL, which will use the hospital beds data for the country provided. If no country is provided then this is 3/100 of hospital beds
`date_ICU_bed_capacity_change`	Calendar dates at which ICU bed capacity changes set in `model_params` change. Defaut = NULL, i.e. no change
`date_vaccine_change`	Date that vaccine doses per day change. Default = NULL.
`baseline_max_vaccine`	Baseline vaccine doses per day. Default = NULL
`max_vaccine`	Time varying maximum vaccine doeses per day. Default = NULL.
`date_vaccine_efficacy_infection_change`	Date that vaccine efficacy against infection changes. Default = NULL.
`baseline_vaccine_efficacy_infection`	Baseline vaccine effacy against infection. Default = NULL
`vaccine_efficacy_infection`	Time varying vaccine efficacy against infection. Default = NULL.
`date_vaccine_efficacy_disease_change`	Date that vaccine efficacy against disease changes. Default = NULL.
`baseline_vaccine_efficacy_disease`	Baseline vaccine efficacy against disease Default = NULL
`vaccine_efficacy_disease`	Time varying vaccine efficacy against infection. Default = NULL.
`Rt_args`	List of arguments to be passed to `evaluate_Rt_pmcmc` for calculating Rt. Current arguments are available in `Rt_args_list`
`burnin`	number of iterations to discard from the start of MCMC run when sampling from the posterior for trajectories
`replicates`	number of trajectories (replicates) to be returned that are being sampled from the posterior probability results produced by `run_mcmc_chain` to select parameter set. For each parmater set sampled, run particle filter with `n_particles` and sample 1 trajectory
`forecast`	Number of days to forecast forward. Default = 0
`required_acceptance_ratio`	Desired MCMC acceptance ratio
`start_adaptation`	Iteration number to begin RM optimisation of scaling factor at
`gibbs_sampling`	Whether or not to use the Gibbs Sampler for start_date
`gibbs_days`	Number of days either side of the start_date parameter to evaluate likelihood at
`...`	Further aguments for the model parameter function. If using the `explicit_model` (default) this will be `parameters_explicit_SEEIR`.

Details

Run a pmcmc sampler with the Squire model setup (i.e. include the various model parameters for the odin model to generate curves)

Value

squire_simulation

output: Trajectories from the sampled pMCMC parameter iterations.
parameters: Model parameters use for squire
model: Squire model used
inputs: Inputs into the squire model for the pMCMC.
pMCMC_results: An mcmc object generated from pmcmc and contains:

inputs: List of inputs
chains: List that include

results: Matrix of accepted parameter samples, rows = iterations as well as log prior, (particle filter estimate of) log likelihood and log posterior
states: Matrix of compartment states
acceptance_rate: MCMC acceptance rate
ess: MCMC chain effective sample size