Core serosolver function running the adaptive Metropolis-within-Gibbs algorithm. It estimates the antibody kinetics parameters and infection histories from the supplied data. The MCMC chains are saved in blocks as CSV files at the location given by filename; the returned object also contains the model settings, diagnostics, and plots when these are requested. See the package guide and case studies for examples.

serosolver(
  par_tab,
  antibody_data,
  demographics = NULL,
  antigenic_map = NULL,
  possible_exposure_times = NULL,
  mcmc_pars = c(),
  mv_proposals = FALSE,
  n_chains = 1,
  parallel = FALSE,
  start_inf_hist = NULL,
  fixed_inf_hists = NULL,
  filename = "test",
  prior_func = NULL,
  prior_version = 2,
  measurement_bias = NULL,
  proposal_ratios = NULL,
  solve_likelihood = TRUE,
  n_alive = NULL,
  random_start_parameters = TRUE,
  start_level = "none",
  data_type = 1,
  verbose = TRUE,
  verbose_dev = FALSE,
  inf_hist_mcmc_summaries = TRUE,
  exponential_waning = FALSE,
  plot_outputs = TRUE,
  ...
)

Arguments

par_tab

The parameter table controlling information such as bounds, initial values etc. See example_par_tab

antibody_data

The data frame of serological measurements to be fitted. It should contain individual, sample_time, biomarker_id, measurement, and birth; biomarker_group, repeat_number, and population_group are optional. If population_group is absent, all individuals are assigned to group 1. See example_antibody_data

demographics

if not NULL, a data frame giving demographic variables for each individual. It must include individual and birth, and must include any variables used to stratify parameters in par_tab. Demographic variables may be fixed for each individual or vary over time. See the demographic stratification and covariate vignette.

antigenic_map

(optional) A data frame of antigenic x and y coordinates. Must have column names: x_coord; y_coord; inf_times. See example_antigenic_map

possible_exposure_times

(optional) this argument gives the vector of times at which individuals can be infected. Defaults to entries in antigenic_map.

mcmc_pars

Named numeric vector with parameters for the MCMC procedure. See details.

mv_proposals

If TRUE, uses a multivariate normal distribution for the proposal distribution. FALSE uses univariate proposals. It is advised to leave this as FALSE, multivariate proposals seems to generally be inefficient for serosolver.

n_chains

Number of MCMC chains to run

parallel

if TRUE, runs multiple chains in parallel using the doParallel package

start_inf_hist

Infection history matrix to start MCMC at. Can be left NULL. See example_inf_hist

fixed_inf_hists

(optional) Data frame with columns individual, time, and value, giving infection states that should be fixed during the MCMC run. See the advanced features vignette.

filename

The file path and prefix for the MCMC output. The chain files are saved with _chain.csv and _infection_histories.csv appended, and the model settings are saved with _serosolver_settings.RData appended. When parallel chains are used, progress messages are written to a matching _log.txt file.

prior_func

User function of prior for model parameters. Should take parameter values only

prior_version

which infection history assumption prior_version to use? See describe_priors for options. Can be 1, 2, 3 or 4

measurement_bias

optional NULL. A data frame mapping each biomarker_id and biomarker_group combination to the rho_index of the measurement-shift parameter that it uses. See the advanced features vignette.

proposal_ratios

optional NULL. Can set the relative sampling weights of the infection state times. Should be an integer vector of length matching the number of infection-history time points. Otherwise, leave as NULL for uniform sampling.

solve_likelihood

if FALSE, returns only the prior and does not solve the likelihood. Use this if you wish to sample directly from the prior

n_alive

if not NULL, uses this as the number alive for the infection history prior, rather than calculating the number alive based on antibody_data

random_start_parameters

if FALSE, uses whatever parameter values were passed in par_tab as the starting positions for the MCMC chain

start_level

either "none" or a data frame giving the starting biomarker level for each individual, biomarker_group, and biomarker_id combination. With "none", starting levels are assumed to be 0. See the advanced features vignette.

data_type

numeric or text value identifying the observation model: 1 or "discrete" for discrete, bounded data; 2 or "continuous" for continuous, bounded data; or 3 or "false_positive" for continuous data with the false-positive observation model. Supply one value per biomarker group, or one value to use for all groups.

verbose

if TRUE, prints progress updates during the run

verbose_dev

if TRUE, prints additional messages regarding step sizes, acceptance rates etc

inf_hist_mcmc_summaries

if TRUE, calculates MCMC summaries of the infection history posterior draws. Set to FALSE to decrease run time, as this is a slow operation.

exponential_waning

Deprecated. If TRUE, assumes exponential waning of antibody levels rather than linear waning, and changes the cross-reactivity model from a linear decline with antigenic distance to an exponential decline. Prefer a fixed exponential_waning row in par_tab, with values = 1 and par_type = 0. See the advanced features vignette.

plot_outputs

if TRUE, calculates diagnostic summaries and returns model-fit, attack-rate, and antibody-model plots. Defaults to TRUE.

...

Other arguments passed to the posterior function.

Value

A list containing the paths to the parameter and infection-history chain files, diagnostic summaries and warnings, the settings used for the fit, fitted antibody predictions, plots when plot_outputs = TRUE, and the loaded MCMC chains.

Details

The mcmc_pars argument is a named vector allowing control over many MCMC options. The key options are:

  • iterations (number of post-adaptation iterations to run)

  • adaptive_iterations (for this many iterations, change proposal step size adaptively every adaptive_frequency iterations)

  • adaptive_frequency (adapt proposal step size every adaptive_frequency iterations)

  • thin (save every n iterations from kinetics parameters samples, advised for long chains, advised to set such that iterations/thin is no more than 1000)

  • thin_inf_hist (save every n iterations from infection history samples, advised to set such that iterations/thin_inf_hist is no more than 1000)

Additional arguments for more precise control, not relevant for the majority of users:

  • proposal_inf_hist_indiv_prop (proportion of individuals resampled at each infection history proposal)

  • save_block (after this many iterations (post thinning), save to disk)

  • target_acceptance_rate_theta (desired acceptance rate for theta parameters)

  • target_acceptance_rate_inf_hist (desired acceptance rate for infection histories)

  • proposal_ratio (ratio of infection history samples to theta samples ie. proposal_ratio = 2 means sample inf hist twice for every theta sample, proposal_ratio = 0.5 means sample theta twice for every inf hist sample)

  • proposal_inf_hist_time_prop (proportion of infection times to resample for each individual at each iteration)

  • proposal_inf_hist_distance (number of infection years/months to move when performing infection history swap step)

  • proposal_inf_hist_adaptive (if 1, performs adaptive infection history proposals. If 0, retains the starting infection history proposal parameters. This should ONLY be turned on for prior version 2)

  • proposal_inf_hist_indiv_swap_ratio (if using gibbs sampling of infection histories, what proportion of proposals should be swap steps, between 0 and 1)

  • proposal_inf_hist_group_swap_ratio (proportion of infection history proposal steps to swap proposal_inf_hist_group_swap_prop of two time periods' contents, between 0 and 1)

  • proposal_inf_hist_group_swap_prop (when swapping contents of two time points, what proportion of individuals should have their contents swapped, between 0 and 1)

  • propose_from_prior (set to 1 to sample directly from the infection history prior, or 0 for independent proposals. Sometimes one version works better than the other, so try switching if you are getting poor infection history convergence)

The parameter and infection-history chains are written to disk during the run, in blocks controlled by save_block. The thin and thin_inf_hist settings control how often parameter and infection-history draws are saved. The settings file records the inputs needed to read the chains back into R with load_mcmc_chains.

Examples

if (FALSE) { # \dontrun{
data(example_antibody_data)
data(example_par_tab)
data(example_antigenic_map)
res <- serosolver(example_par_tab, example_antibody_data, example_antigenic_map)
} # }