| Title: | Plague Transmission Models for Epidemiological Research |
|---|---|
| Description: | A toolkit for plague transmission modeling using stochastic compartmental models. Implements carcass-based transmission dynamics (Didelot et al. 2017) with spatial structure and multi-host dynamics for epidemiological research. |
| Authors: | Nick Gauthier [aut, cre] (ORCID: <https://orcid.org/0000-0002-2225-5827>), Nich Martin [aut] (ORCID: <https://orcid.org/0000-0001-9539-250X>) |
| Maintainer: | Nick Gauthier <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2026-06-03 22:13:24 UTC |
| Source: | https://github.com/flmnh-ai/yersinia |
Mirrors the canonical list in load_named_scenario(). v1 hardcodes the
five built-in scenarios; future versions will scan
system.file("scenarios", package = "yersinia") so user-bundled YAMLs
show up automatically.
available_scenarios()available_scenarios()
Character vector of scenario names.
Vectorized one-step update of a 3-state CA (0 = susceptible, 1 =
epidemic, 2 = endemic) following K&G (2000) Table 3 semantics:
ca_step( state, sus, A, lookups, d_per_tick, P_recovery_sus = 20, grow_per_tick = 12 )ca_step( state, sus, A, lookups, d_per_tick, P_recovery_sus = 20, grow_per_tick = 12 )
state |
Integer vector of length |
sus |
Numeric vector of length |
A |
Sparse adjacency matrix from |
lookups |
List of three functions ( |
d_per_tick |
Probability per tick that an isolated |
P_recovery_sus |
Susceptibility (0-100) that recovered |
grow_per_tick |
Increment in |
S → E or S → P if any neighbour is in E or P, with per-target
probability 1 - (1 - T_E)^n_E × (1 - T_P)^n_P and outcome P with
conditional probability Q (else E).
E → S_0 (deterministic — epidemics resolve in one tick).
P → S_w with probability d_per_tick, sus reset to P_recovery_sus.
Susceptibility regrows in S cells by grow_per_tick (capped at 100).
All transitions are computed from the start-of-tick state; no within-tick asynchronous updates.
A list with the updated state and sus vectors.
Calculate equilibrium values for deterministic models
calculate_equilibrium(params, model_type = "rats_only")calculate_equilibrium(params, model_type = "rats_only")
params |
plague_parameters object or list |
model_type |
Type of model |
List of equilibrium values
This function is for reference/educational purposes. For operational models, use the stochastic simulations in run_plague_model().
Calculate force of infection over time
calculate_force_of_infection(results)calculate_force_of_infection(results)
results |
plague_results object |
Tibble with force of infection by time
Calculate outbreak summary statistics
calculate_outbreak_metrics(results, compartment = "I", threshold = 1)calculate_outbreak_metrics(results, compartment = "I", threshold = 1)
results |
plague_results object |
compartment |
Compartment to analyze (default "I") |
threshold |
Minimum value to consider as "outbreak" (default 1) |
Data frame with outbreak statistics
Calculate basic reproduction number (R0) for different model types
calculate_R0(params, model_type = "rats_only", K_r = 2500, K_h = 5000)calculate_R0(params, model_type = "rats_only", K_r = 2500, K_h = 5000)
params |
plague_parameters object or list |
model_type |
Type of model ("rats_only", "with_humans", etc.) |
K_r |
Rat carrying capacity (defaults to 2500 if not in params) |
K_h |
Human carrying capacity (defaults to 5000 if not in params, only used for "with_humans") |
Numeric R0 value or list of R0 values
For multi-outbreak fits, dust2's grouped likelihood requires every group
to share the same time grid; shorter outbreaks are NA-padded out to the
longest. The returned tibble has columns group (character — dust2
rejects factor group), time, deaths.
cohort_data(cohort_ids, data = NULL)cohort_data(cohort_ids, data = NULL)
cohort_ids |
Character vector of |
data |
Long outbreaks tibble (defaults to bundled outbreaks). |
Tibble with group, time, deaths.
dust2's grouped likelihood evaluates a single obs_period; mixing daily
and weekly outbreaks in one fit isn't supported here (London 1563 is
weekly; everything else is daily). Errors when the cohort mixes cadences.
cohort_obs_period(cohort_ids, data = NULL)cohort_obs_period(cohort_ids, data = NULL)
cohort_ids |
Character vector of |
data |
Long outbreaks tibble (defaults to bundled outbreaks). |
Single integer.
Per-outbreak population values (K_h / K_r when fixed).
cohort_population(cohort_ids, data = NULL)cohort_population(cohort_ids, data = NULL)
cohort_ids |
Character vector of |
data |
Long outbreaks tibble. |
Named numeric, one entry per outbreak.
Selection state is lab_session$cohort_ids (character vector of
outbreak_ids). The UI is rebuilt whenever that vector changes.
cohort_server(id, lab_session, data = NULL)cohort_server(id, lab_session, data = NULL)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
data |
Long outbreaks tibble. Defaults to the bundled outbreaks. |
Invisibly, the moduleServer result.
Cohort module — UI.
cohort_ui(id)cohort_ui(id)
id |
Module namespace id. |
A bslib::card() with the cohort chips + picker grid.
Excludes parameters that aren't user-fittable in practice:
configurable_param_names()configurable_param_names()
tau (timestep — system, not biological)
seasonal (vector — handled separately as climate input)
obs_period (integer flag — derived from the cohort data)
iota (resistance fecundity multiplier — always-fixed per CLAUDE.md)
The returned set is ~23 scalar parameters: rates, capacities, initial
counts, and the likelihood dispersion kappa.
Character vector of parameter names.
Returns a Brière-shaped multiplier capturing the combined temperature effects on flea survival and Y. pestis biofilm blockage formation. Multiply your fitted baseline delta_R rate by this output to obtain delta_R(T) for use in the model:
delta_R_multiplier( temperature, T_min = 5, T_max = 37, q = 1.5, T_ref = 17.8, cap = 1000 )delta_R_multiplier( temperature, T_min = 5, T_max = 37, q = 1.5, T_ref = 17.8, cap = 1000 )
temperature |
Numeric vector of temperatures (degrees Celsius) |
T_min |
Lower thermal bound (default 5 C) |
T_max |
Upper thermal bound (default 37 C) |
q |
Warm-side asymmetry exponent (default 1.5) |
T_ref |
Reference temperature where multiplier = 1 (default 17.8 C) |
cap |
Maximum multiplier value, applied beyond thermal bounds (default 1000) to keep the likelihood numerically well-behaved |
delta_R(T) = delta * delta_R_multiplier(T)
Functional form: L(T) = (T - T_min) * (T_max - T)^q for T_min < T < T_max multiplier(T) = L(T_ref) / L(T)
Defaults place the optimum (longest infectious period) at 17.8 C with hard bounds at 5 C (cold cliff: Y. pestis blockage formation fails) and 37 C (warm cliff: flea thermal death + biofilm dissolution). These are anchored to: Krauer et al. 2021 (epidemic peak ~17 C), Hinnebusch lab biofilm work (cold cliff), and Mellanby 1932 (thermal death point of X. cheopis).
Default T_ref = 17.8 means the fitted delta represents delta_R at the
thermal optimum, i.e. the minimum loss rate (longest carcass
infectiousness). All multipliers are then >= 1. To match Didelot et al.
2017's calibration anchor instead, set T_ref = 14 — delta then
represents delta_R at Cairo winter mean temperature, comparable to
their posterior median of 0.267/day.
Numeric vector of multipliers, one per input temperature
Diagnostics strip module — server.
diag_strip_server(id, lab_session)diag_strip_server(id, lab_session)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
Invisibly, the moduleServer result.
Diagnostics strip module — UI.
diag_strip_ui(id)diag_strip_ui(id)
id |
Module namespace id. |
A shiny::div() for the strip (initially empty).
For each parameter with finite [low, high] bounds in bounds, computes
the fraction of draws within edge (a fraction of the bound width) of
either bound. Flags any parameter where that fraction exceeds threshold.
Bound-piling indicates the prior is constraining the posterior — either
the prior is too tight, or the model wants a parameter the prior forbids.
diagnose_bound_piling(samples, bounds, threshold = 0.2, edge = 0.05)diagnose_bound_piling(samples, bounds, threshold = 0.2, edge = 0.05)
samples |
A |
bounds |
Named list keyed by parameter name; each entry is a length-2
numeric |
threshold |
Fraction of draws within |
edge |
Fraction of bound width counted as "near the edge". Default 0.05 (5%). |
Parameters with non-finite bounds (e.g. Exponential(rate=0.1) with
high = Inf) are silently skipped — there's no upper edge to pile against.
A diagnostic record or NULL.
Uses posterior::rhat() on each parameter; flags the fit if any parameter's
R-hat exceeds threshold. The MVP threshold of 1.5 is intentionally loose
(the conventional convergence cutoff is 1.05–1.1) — we want to flag pilots
that are obviously broken, not over-warn on borderline cases.
diagnose_chain_stuck(samples, threshold = 1.5)diagnose_chain_stuck(samples, threshold = 1.5)
samples |
A |
threshold |
R-hat above which the fit is flagged. Default 1.5. |
A diagnostic record (see file header) or NULL if all R-hats are
below threshold.
kappa).The plague-fit likelihood is deaths ~ NegBinomial(size = kappa, mu = ...).
Small kappa means highly overdispersed observations — the data is so
noisy that almost any trajectory is plausible, so the posterior over
transmission parameters becomes uninformative. The MVP threshold of 5 is
a rule of thumb from the Cairo / Barcelona pilots in this conversation:
kappa < 5 typically corresponded to broken or under-specified fits.
diagnose_low_kappa(samples, par = "kappa", threshold = 5)diagnose_low_kappa(samples, par = "kappa", threshold = 5)
samples |
A |
par |
Name of the dispersion parameter. Default |
threshold |
Median |
Uses the posterior median of kappa across all chains and iterations.
A diagnostic record or NULL (also NULL if par isn't in the
sampled parameters — kappa may be pinned in some configurations).
Renders cohort death curves and (when a fit exists) the posterior
predictive trajectory fan. Recomputes the fan lazily on fit_state
changes.
hero_server(id, lab_session, n_draws = 80L)hero_server(id, lab_session, n_draws = 80L)
id |
Module namespace id (must match |
lab_session |
A LabSession instance. |
n_draws |
Number of posterior draws to render. Default 80. |
Invisibly, the moduleServer result.
Hero plot module — UI.
hero_ui(id)hero_ui(id)
id |
Module namespace id. |
A bslib::card() with the hero plot.
Subsamples n_draws posterior draws (evenly spaced by posterior::thin_draws),
runs the deterministic plague humans model forward for each, and returns a
long-format tibble of per-draw predicted observation series. Used by the
hero plot to render the posterior fan over cohort data.
lab_fit_forward_sim(setup, samples, n_draws = 100, t_grid = NULL)lab_fit_forward_sim(setup, samples, n_draws = 100, t_grid = NULL)
setup |
Output of |
samples |
A |
n_draws |
Target number of posterior draws to simulate. Default 100. |
t_grid |
Optional integer vector of simulation times. Defaults to
|
Handles both the single-outbreak (flat packer) and multi-outbreak (grouped packer) cases.
Tibble with columns draw, group, time, mu (predicted
observation expectation per posterior draw).
lab_fit_assemble().Synchronous: blocks until all chains complete. v1 uses
monty::monty_runner_serial() (sequential chains in the calling process)
to keep the runtime model simple; switch to monty_runner_callr in v2
for parallel chains.
lab_fit_run( setup, n_chains = 4L, n_iter = 30000L, burnin_frac = 0.5, initial = NULL )lab_fit_run( setup, n_chains = 4L, n_iter = 30000L, burnin_frac = 0.5, initial = NULL )
setup |
Output of |
n_chains |
Number of chains. Default 4. |
n_iter |
Iterations per chain. Default 30000. Single-outbreak fits
converge by ~10k; multi-outbreak hierarchical fits (4+ outbreaks)
typically need 30–50k for the adaptive sampler to settle the
per-group |
burnin_frac |
Fraction of iterations to discard from each chain as warmup. Default 0.5 (drop first half). Set to 0 to keep everything. |
initial |
Optional matrix of initial parameter vectors
( |
Discards the first burnin_frac of iterations from each chain post-hoc.
The adaptive sampler tunes its proposal VCV continuously during sampling,
so early iterations are biased relative to the converged regime —
discarding the first half is the convention for adaptive RWM (see e.g.
Stan / NUTS warmup defaults). Returns the trimmed samples.
A monty_samples object with the warmup iterations removed.
After a pilot completes, use the adaptive sampler's tuned proposal VCV to warm-start a shorter stochastic chain that respects the actual particle-filter likelihood (the pilot's deterministic unfilter is an expected-value approximation). Returns the production samples; the caller is responsible for keeping the pilot around for diagnostics.
lab_fit_run_production( pilot_setup, pilot_samples, n_chains = 4L, n_iter = 5000L, burnin_frac = 0.2, n_particles = 500L, n_threads = 2L, vcv_inflation = NULL )lab_fit_run_production( pilot_setup, pilot_samples, n_chains = 4L, n_iter = 5000L, burnin_frac = 0.2, n_particles = 500L, n_threads = 2L, vcv_inflation = NULL )
pilot_setup |
The setup list returned by |
pilot_samples |
The samples returned by |
n_chains |
Number of chains. Default 4. |
n_iter |
Iterations per chain. Default 5000 (much shorter than pilot since the proposal is already tuned). |
burnin_frac |
Fraction to discard. Default 0.2 (chains start from a converged pilot region, so minimal warmup needed). |
n_particles |
Particle filter size. Default 500. |
n_threads |
OpenMP threads. Default 2. |
vcv_inflation |
Factor to scale the tuned pilot VCV by. Default
|
Mirrors the two-stage workflow in
vignettes/monty-barcelona-hierarchical.qmd and friends.
A monty_samples object with the warmup iterations removed.
Inverse of LabSession$to_list(). Unknown keys in state are ignored so
older snapshots remain loadable as the schema grows; missing keys fall
back to the same defaults LabSession's constructor uses.
lab_session_from_list(state)lab_session_from_list(state)
state |
Named list as returned by |
A new LabSession instance.
LabSession R6 class — Virtual Lab session state.
LabSession R6 class — Virtual Lab session state.
Holds the reactive state of an interactive plague-fitting session. Used by
lab_app() and the Shiny modules. State is exposed as active bindings
backed by shiny::reactiveVal() so reads from inside reactive contexts
subscribe automatically; outside Shiny, reads return the current value.
Current fields: cohort_ids (character vector of selected outbreak_ids
from the outbreaks dataset), model_config (list with scenario,
shared, local — see model_config_default()). Future versions will
add priors, fit_state, sandbox_draw as the corresponding modules
come online.
Use lab_session_from_list() to reconstruct from a saved state list.
cohort_idsCharacter vector of selected outbreak_ids.
Reading from inside a reactive context subscribes; assigning
triggers downstream observers.
model_configModel-config list (scenario, shared, local)
— see model_config_default() for the schema.
priorsNamed list of priors keyed by parameter name; each
entry is list(family, params). See priors_default().
fit_stateNamed list summarising the most recent fit:
status ("idle" | "running" | "complete" | "error"),
samples (a monty_samples object or NULL),
n_chains, n_iter, duration (seconds), error (message),
started_at, completed_at.
new()
Construct a new LabSession.
LabSession$new( cohort_ids = character(0), model_config = NULL, priors = NULL, fit_state = NULL )
cohort_idsInitial cohort (character vector of outbreak_ids).
Default empty.
model_configInitial model config; defaults to
model_config_default().
priorsNamed list of priors. Defaults to priors_default()
over the fitted parameters implied by model_config.
fit_stateInitial fit-state struct; defaults to an idle state.
to_list()
Snapshot the session state as a plain list (no reactive
infrastructure), suitable for saveRDS().
LabSession$to_list()
apply_list()
Mutate this LabSession in place from a saved-state list
(e.g. one produced by to_list()). All reactiveVals fire, so any
subscribed Shiny observers update. Unknown keys are ignored so
older snapshots remain loadable.
LabSession$apply_list(state)
stateNamed list as returned by to_list().
clone()
The objects of this class are cloneable with this method.
LabSession$clone(deep = FALSE)
deepWhether to make a deep clone.
Load named scenario parameters
load_named_scenario(name)load_named_scenario(name)
name |
Name of scenario |
List of model parameters
Load plague model scenario parameters
load_scenario(scenario = "defaults", ...)load_scenario(scenario = "defaults", ...)
scenario |
Name of scenario or path to YAML file or list of parameters |
... |
Additional model parameters to override |
List of model parameters
Constructs a hex-tessellated grid covering a side × side bounding box
(one hex cell per unit area, pointy-top by default), with 6-neighbour
adjacency. Requires the sf package. Returns the same lattice interface
as make_square_grid() plus an sf field that lets you visualize results
with ggplot2::geom_sf().
make_hex_lattice(side)make_hex_lattice(side)
side |
Numeric. Side length of the square bounding box. The number
of cells produced is approximately |
A list with neighbors, A, npop, and sf (an sf data frame
with cell_id and geometry columns).
make_square_grid() for the square variant.
From a tidy data frame of empirical T_i / Q_i estimates (as produced
by the metapop-cellular-automaton vignette), build linear interpolation
functions for T_E(S_frac), T_P(S_frac), and Q(S_frac) at a chosen
coupling level mu_r. The returned lookup list is the form
ca_step() expects.
make_lookups(TQ, mu_r_use)make_lookups(TQ, mu_r_use)
TQ |
Data frame with columns |
mu_r_use |
Coupling level at which to extract lookups (must match a
value in |
Pads T_i only with (S_frac = 0, T_i = 0) (biologically required —
fully-resistant cells cannot trigger) and uses approxfun(rule = 2) to
extrapolate T_i past the empirical S_frac range as the boundary value.
For Q_i no synthetic boundary points are added, since the empirical
Q curve often drops at high S_frac (fully-susceptible patches burn out
fast) and a Q(1) = 1 pad would inflate the lookup exactly where the
data says it should fall. Rare-trigger NAs in Q_i (n_trig < 5) are
dropped before interpolation.
A list of function(S_frac) interpolators: T_E, T_P, Q.
Constructs a n_row × n_col regular grid with rook (4-neighbour)
adjacency, returning a list with the neighbour list, sparse adjacency
matrix, and per-cell (row, col) coordinates. Boundary cells have 2-3
neighbours; interior cells have 4. The neighbour ordering is row-major
(cell index = (row - 1) * n_col + col).
make_square_grid(n_row, n_col)make_square_grid(n_row, n_col)
n_row, n_col
|
Integer dimensions of the lattice. |
A list with elements:
neighbors — list of integer vectors, one per cell
A — sparse Matrix::sparseMatrix adjacency, dims (npop, npop)
n_row, n_col, npop
cell_xy — data frame with row and col columns
make_hex_lattice() for the hex variant; run_ca() consumes
the returned list directly.
Model accordion module — server.
model_accordion_server(id, lab_session)model_accordion_server(id, lab_session)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
Invisibly, the moduleServer result.
Model accordion module — UI.
model_accordion_ui(id)model_accordion_ui(id)
id |
Module namespace id. |
A bslib::accordion() containing the model config controls.
Mirrors the Barcelona / Cairo vignette fitting workflows: biology +
likelihood-overdispersion are shared across the cohort; transmission to
humans and ignition + record-quality knobs are local (per outbreak).
Populations (K_h, K_r) are not in the fitted set — when omitted,
lab_fit_assemble() pins them per-group at each outbreak's known
population from the outbreaks dataset. Same applies to g_h — the
vignettes pin it at the scenario value because it's only weakly
identifiable from death counts alone.
model_config_default(scenario = "defaults")model_config_default(scenario = "defaults")
scenario |
Base scenario name. Default |
rho (flea redistribution exponent) is also intentionally not fitted
by default: it's confounded with beta_r via the rat-side R₀ formula
R₀ ∝ beta_r · (1 − e^(−rho)) / delta_R, so jointly fitting both
creates a ridge in the posterior that wrecks rhat. Pin rho at the
scenario value and let beta_r carry the rat-side identification.
To fit K_h/K_r, g_h, or rho, add them via the Model tab.
A model-config list.
Loads the named scenario, then for every configurable parameter reports
the scenario value, the chosen scope (fixed | shared | local), and the
value the model would actually see (scenario_value for fixed, NA
for shared/local since the sampler will set those).
model_config_resolve(config)model_config_resolve(config)
config |
A model-config list (see file header). |
Parameters that appear in both shared and local are reported as
local — the multi-cohort case usually wants per-outbreak resolution
when the user has flagged conflict.
A tibble with columns parameter, scope, scenario_value,
resolved_value, sorted alphabetically.
Create a plague_results object
new_plague_results(data, model_type, params, run_info = list())new_plague_results(data, model_type, params, run_info = list())
data |
Tibble with simulation results |
model_type |
Type of model used |
params |
Parameters used in simulation |
run_info |
Runtime information |
plague_results object
Aggregates the long-format outbreaks dataset to one row per outbreak, with the metadata + summary stats the cohort picker displays (location, year, population, duration, total deaths, peak day, observation period).
outbreak_summary(data = NULL)outbreak_summary(data = NULL)
data |
Long-format outbreaks tibble. Defaults to the package's bundled outbreaks dataset. |
A tibble with one row per outbreak_id, sorted by year.
outbreak_summary()outbreak_summary()
Historical plague mortality data from multiple European and Middle Eastern cities spanning from 1348 to 1835, compiled by Dean et al. (2018).
outbreaksoutbreaks
A tibble with 1,451 rows and 7 variables:
Unique identifier for each outbreak (location_year)
City name
Year of the outbreak
Total population of the city
Calendar day of the observation, 1-indexed from outbreak start. For weekly outbreaks (London 1563) this is the day at the end of the reporting week (7, 14, ...).
Death count for the reporting window ending at day;
daily for most outbreaks, weekly for London 1563 – see caveats
below re: plague-specific vs all-cause across sources
Length in days of the reporting window each
deaths value covers (1 for daily, 7 for London 1563). Pass
this to plague_fit_setup via the obs_period
argument so the model's D_h accumulator aggregates over a
matching window.
The dataset includes plague outbreaks from:
Barcelona 1490 (daily data, 182 days)
Malta 1813 (daily data, 209 days)
Florence 1400 (daily data, 180 days)
Cairo 1835 (daily data, 366 days)
Eyam 1665 (daily data, 145 days)
Prague 1713 (daily data, 198 days)
London 1563 (weekly data, 33 weeks; obs_period = 7)
Givry 1348 (daily data with some missing values, 138 days)
Some datasets contain missing values (NA) where historical records were
incomplete. London 1563 carries weekly Bills of Mortality counts – the
obs_period = 7 flag signals that to fitting code, and day
steps in increments of 7.
Cause-of-death provenance. The source documents differ in whether
they record plague-specific deaths or all burials during the outbreak. This
matters when fitting to the model's D_h (plague deaths only). Summary:
Plague-specific (cause distinguished in source): Barcelona 1490 (cerca de morts), London 1563 (Bills of Mortality), Malta 1813, Prague 1713, Cairo 1835
All-cause burial records, but with negligible baseline vs. epidemic peak: Givry 1348 (parish register, ~1,200 pop), Eyam 1665 (parish register, ~700 pop). Fittable as plague-specific in practice.
Ambiguous: Florence 1400 (Libri dei Morti, inconsistent cause-of-death annotation)
See data-raw/outbreaks.R for per-outbreak primary-source citations.
Users can filter to specific outbreaks using the outbreak_id or location/year:
# Get Barcelona data
barcelona <- plague_outbreaks |> filter(outbreak_id == "Barcelona_1490")
# Get all 15th century outbreaks
fifteenth_century <- plague_outbreaks |> filter(year >= 1400 & year < 1500)
Dean, K.R., Krauer, F., Walløe, L., Lingjærde, O.C., Bramanti, B., Stenseth, N.C. and Schmid, B.V., 2018. Human ectoparasites and the spread of plague in Europe during the Second Pandemic. Proceedings of the National Academy of Sciences, 115(6), pp.1304-1309.
# View structure of the data str(outbreaks) # Summary statistics by outbreak outbreaks |> group_by(outbreak_id, location, year) |> summarise( duration_days = max(day), total_deaths = sum(deaths, na.rm = TRUE), peak_deaths = max(deaths, na.rm = TRUE), population = first(population) )# View structure of the data str(outbreaks) # Summary statistics by outbreak outbreaks |> group_by(outbreak_id, location, year) |> summarise( duration_days = max(day), total_deaths = sum(deaths, na.rm = TRUE), peak_deaths = max(deaths, na.rm = TRUE), population = first(population) )
Build a particle filter for the humans plague model.
plague_fit_filter(data, n_particles = 2000, n_threads = 2)plague_fit_filter(data, n_particles = 2000, n_threads = 2)
data |
Data frame with columns |
n_particles |
Number of particles in the filter. |
n_threads |
Number of OpenMP threads for the filter (default 2). |
A dust2 particle filter ready for
dust2::dust_likelihood_run() or dust2::dust_likelihood_monty().
These are the parameters swept by the monty sampler in plague_fit_setup().
If you change this set, you also need to provide a matching prior and
vcv — the defaults returned by plague_fit_prior() and plague_fit_vcv()
are built for exactly this ordering.
plague_fit_fitted_names()plague_fit_fitted_names()
The model also supports p_obs (observation/ascertainment probability,
defaults to 1 = perfect reporting) which is deliberately NOT fitted by
default. To opt in, pass e.g.
fitted_names = c(plague_fit_fitted_names(), "p_obs") to
plague_fit_setup() with a matching prior
(p_obs ~ Uniform(0, 1) follows the 2009-flu example in the odin-monty
book) and an extended VCV.
Character vector of parameter names.
Loads a scenario via load_scenario(), drops any parameter that the
sampler will vary, and filters to just the keys the odin2 humans model
accepts. The returned list is suitable for the fixed argument of
monty::monty_packer() or as pars for dust2::dust_system_create().
plague_fit_fixed_pars( scenario = "historical", fitted_names = plague_fit_fitted_names() )plague_fit_fixed_pars( scenario = "historical", fitted_names = plague_fit_fitted_names() )
scenario |
Scenario name, path to YAML, or parameter list. Defaults
to |
fitted_names |
Names of parameters that the sampler will set; these
are stripped from the fixed set. Defaults to |
Plain named list of fixed parameters.
Paired with plague_fit_fitted_names() — order and variables must match.
plague_fit_prior()plague_fit_prior()
A monty model returned by monty::monty_dsl().
monty::monty_sample() on a plague fit.Wires the filter, packer (with fixed parameters), likelihood, prior, and
sampler together. Runs the filter once at the fixed defaults as a smoke
test — if that throws, the parameters and data don't agree and you want
the error now rather than inside monty_sample().
plague_fit_setup( data, scenario = "historical", fitted_names = plague_fit_fitted_names(), n_particles = 2000, n_threads = 2, prior = NULL, vcv = NULL, obs_period = 1L, smoke_test = TRUE )plague_fit_setup( data, scenario = "historical", fitted_names = plague_fit_fitted_names(), n_particles = 2000, n_threads = 2, prior = NULL, vcv = NULL, obs_period = 1L, smoke_test = TRUE )
data |
Data frame with columns |
scenario |
Fixed-parameter scenario (name, path, or list). Default
|
fitted_names |
Parameters for the sampler to vary. Default
|
n_particles |
Particle filter size. Default 2000. |
n_threads |
OpenMP threads for the filter. Default 2. |
prior |
monty model; default |
vcv |
Sampler VCV; default |
obs_period |
Integer length (in days) of the observation window over
which |
smoke_test |
If |
Named list: posterior, sampler, filter, packer,
fixed_pars, fitted_names. The caller chooses runner, n_steps,
n_chains, burnin when invoking monty_sample().
Diagonal entries are ordered to match plague_fit_fitted_names().
Replace with an empirical covariance after a pilot run for better mixing
(the monty vignette walks through that tuning step).
plague_fit_vcv()plague_fit_vcv()
Numeric matrix.
An odin2/dust2 system generator compiled at package build time from
inst/odin/plague_stochastic_humans.R. Pass it directly to
dust2::dust_system_create() or dust2::dust_filter_create().
An odin2/dust2 system generator compiled at package build time from
inst/odin/plague_stochastic_metapop.R. Same Didelot carcass-based
transmission core as plague_stochastic_humans, with rat S/I/R compartments
arrayed over npop patches and a row-stochastic contact_r[npop, npop]
matrix routing migrants. Carcasses Q are sessile; humans live per-patch
and do not migrate (v1). Pass directly to dust2::dust_system_create().
Create multi-panel comparison plot
plot_comparison(results_list, compartment = "I", population = 1)plot_comparison(results_list, compartment = "I", population = 1)
results_list |
Named list of plague_results objects |
compartment |
Compartment to compare (default "I") |
population |
Which population to plot (default 1) |
ggplot2 object
Plot compartment dynamics with uncertainty bands
plot_dynamics( results, compartments = NULL, population = 1, show_uncertainty = TRUE, log_scale = FALSE )plot_dynamics( results, compartments = NULL, population = 1, show_uncertainty = TRUE, log_scale = FALSE )
results |
plague_results object |
compartments |
Vector of compartments to plot (NULL for all) |
population |
Which population to plot (for spatial models, default 1) |
show_uncertainty |
Show uncertainty bands for stochastic models (default TRUE) |
log_scale |
Use log scale for y-axis (default FALSE) |
ggplot2 object
Plot phase portrait (S vs I)
plot_phase_portrait( results, compartments = c("S", "I"), population = 1, replicate = 1, ... )plot_phase_portrait( results, compartments = c("S", "I"), population = 1, replicate = 1, ... )
results |
plague_results object |
compartments |
Vector of two compartments to plot (default c("S", "I")) |
population |
Which population to plot (for spatial models, default 1) |
replicate |
Which replicate to plot (for stochastic models, default 1) |
... |
Additional ggplot2 arguments |
ggplot2 object
Plot parameter sensitivity results
plot_sensitivity(sensitivity_results, compartment = "I", metric = "peak")plot_sensitivity(sensitivity_results, compartment = "I", metric = "peak")
sensitivity_results |
Output from run_sensitivity_analysis |
compartment |
Compartment to plot (default "I") |
metric |
Summary metric to plot ("peak", "final", "auc") |
ggplot2 object
Plot method for plague_results
## S3 method for class 'plague_results' plot(x, compartments = NULL)## S3 method for class 'plague_results' plot(x, compartments = NULL)
x |
plague_results object |
compartments |
Vector of compartments to plot (NULL for all) |
Print method for outbreak metrics
## S3 method for class 'plague_outbreak_metrics' print(x, ...)## S3 method for class 'plague_outbreak_metrics' print(x, ...)
x |
plague_outbreak_metrics object |
... |
Additional arguments (ignored) |
Print method for plague_results
## S3 method for class 'plague_results' print(x, ...)## S3 method for class 'plague_results' print(x, ...)
x |
plague_results object |
... |
Additional arguments passed to tibble print |
Print method for plague_parameters
## S3 method for class 'scenario_parameters' print(x, ...)## S3 method for class 'scenario_parameters' print(x, ...)
x |
plague_parameters object |
... |
Additional arguments (ignored) |
Returns a length-2 numeric c(low, high). Only families with a finite
interval support produce both ends finite — Uniform and Beta. For
families with infinite or one-sided support (Normal, LogNormal,
Gamma, Exponential), the unbounded side is Inf/-Inf. The
diagnose_bound_piling() detector skips parameters with non-finite
bounds, so this function safely "no-ops" the diagnostic for unbounded
priors.
prior_bounds(prior)prior_bounds(prior)
prior |
A prior list ( |
Length-2 numeric c(low, high).
Hand-picked defaults seeded from the literature and from the pilots in the
existing monty vignettes. Returns NULL when there is no canonical default
for parameter — the caller should fall back to Uniform(0, 1) or prompt
the user.
prior_default(parameter)prior_default(parameter)
parameter |
Name of a model parameter. |
A prior list (family, params) or NULL.
Returns a tibble with columns x and density, suitable for ggplot.
The grid spans the family's prior_families() range(params) with n
points; densities outside the support of compactly-supported families
(Uniform, Beta) are zero by definition of the underlying R function.
prior_density(prior, n = 200)prior_density(prior, n = 200)
prior |
A prior list ( |
n |
Number of grid points. Default 200. |
A tibble with columns x and density.
Each entry is a list with:
params: named list of parameter specs (default, min, max, label),
in the order the family expects them when called via the monty DSL
density: function (x, p) returning the PDF at x, where p is a
named list of parameter values
range: function (p) returning a length-2 numeric c(low, high) —
sensible plot range for the family at parameter values p
prior_families()prior_families()
v1 supports six families: Uniform, Normal, LogNormal, Gamma,
Exponential, Beta. Names match the monty DSL exactly.
Named list of family definitions.
Produces an unevaluated expression of the form name ~ Family(args),
where the args are taken from prior$params in the order the family
declares them. Suitable for splicing into monty::monty_dsl().
prior_to_dsl(name, prior)prior_to_dsl(name, prior)
name |
Parameter name (will become the LHS of |
prior |
A prior list ( |
A length-1 unevaluated expression.
Priors accordion module — server.
priors_accordion_server(id, lab_session)priors_accordion_server(id, lab_session)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
Invisibly, the moduleServer result.
Priors accordion module — UI.
priors_accordion_ui(id)priors_accordion_ui(id)
id |
Module namespace id. |
A bslib::card() whose body is the dynamic priors accordion.
For each parameter in fitted_params, returns its prior_default() or
falls back to Uniform(0, 1) if no canonical default exists.
priors_default(fitted_params)priors_default(fitted_params)
fitted_params |
Character vector of parameter names. |
Named list of priors keyed by parameter name.
Maps each packer_names entry back to its base parameter (stripping any
<group> decoration), looks up that parameter's prior, and returns the
bounds via prior_bounds(). Suitable for passing to
diagnose_bound_piling().
priors_to_bounds(priors, packer_names)priors_to_bounds(priors, packer_names)
priors |
Named list of priors. |
packer_names |
Character vector of (possibly decorated) parameter
names from |
Named list of c(low, high) keyed by packer_names.
Convert a full priors list to a list of monty-DSL expressions.
priors_to_dsl(priors)priors_to_dsl(priors)
priors |
Named list of priors (see file header). |
List of unevaluated expressions in the same order as priors.
n_seeds seed cells for n_ticks
Initializes a fresh lattice (all S with full susceptibility), seeds
n_seeds cells in state E, and advances ca_step() for n_ticks.
Returns the per-tick state matrix; the sus history is not retained
(regenerable from state and the recovery rules if needed).
run_ca(lattice, lookups, d_per_tick, n_ticks, n_seeds = 5, seed = NULL, ...)run_ca(lattice, lookups, d_per_tick, n_ticks, n_seeds = 5, seed = NULL, ...)
lattice |
Lattice list from |
lookups |
Lookup list from |
d_per_tick |
|
n_ticks |
Integer number of ticks to advance. |
n_seeds |
Integer number of cells to seed in state |
seed |
Optional integer for |
... |
Additional arguments forwarded to |
Integer matrix of shape (n_ticks + 1, npop) with row t + 1
the state at tick t.
Bundles diagnose_chain_stuck(), diagnose_bound_piling(), and
diagnose_low_kappa(). Returns a list of records (length 0 to 3); pass
bounds = NULL to skip bound-piling.
run_diagnostics( samples, bounds = NULL, rhat_threshold = 1.5, bound_threshold = 0.2, bound_edge = 0.05, kappa_par = "kappa", kappa_threshold = 5 )run_diagnostics( samples, bounds = NULL, rhat_threshold = 1.5, bound_threshold = 0.2, bound_edge = 0.05, kappa_par = "kappa", kappa_threshold = 5 )
samples |
A |
bounds |
Optional named list of |
rhat_threshold |
Passed to |
bound_threshold |
Passed to |
bound_edge |
Passed to |
kappa_par |
Passed to |
kappa_threshold |
Passed to |
List of diagnostic records (possibly empty).
Run human stochastic model (dust2/odin2 backend)
run_human_stochastic_model(params, timesteps, n_particles, n_threads)run_human_stochastic_model(params, timesteps, n_particles, n_threads)
params |
List of parameters |
timesteps |
Vector of timesteps |
n_particles |
Number of particles |
n_threads |
Number of threads |
Tidy tibble with results
Forward-simulates the metapopulation plague model
(inst/odin/plague_stochastic_metapop.R). Same Didelot carcass-based
transmission core as run_plague_model(), with rats dispersing between
patches at rate mu_r (via contact_r) and humans dispersing at rate
mu_h (via contact_h). For each species, susceptible, infected, and
recovered individuals migrate at the same rate; carcasses Q are sessile.
Human-side: I_h represents both incubating and symptomatic infections
(the model has no E_h compartment), and small mu_h is the parsimonious
way to express that humans are sedentary on outbreak timescales.
run_plague_metapop_model( scenario = "defaults", npop, contact_r, mu_r, contact_h = NULL, mu_h = 0, K_r, K_h, I_ini, R_ini = NULL, I_h_ini = NULL, R_h_ini = NULL, years = 1, timestep = c("daily", "weekly"), n_particles = 100, n_threads = 1, obs_period = 1L, seasonal = NULL, ... )run_plague_metapop_model( scenario = "defaults", npop, contact_r, mu_r, contact_h = NULL, mu_h = 0, K_r, K_h, I_ini, R_ini = NULL, I_h_ini = NULL, R_h_ini = NULL, years = 1, timestep = c("daily", "weekly"), n_particles = 100, n_threads = 1, obs_period = 1L, seasonal = NULL, ... )
scenario |
Parameter set name, file path, or list of parameters. |
npop |
Integer number of patches. |
contact_r |
Row-stochastic |
mu_r |
Rat migration rate per day. Same value applies to S, I, R rats. |
contact_h |
Row-stochastic |
mu_h |
Human migration rate per day (default 0 – sessile humans). Same value applies to S_h, I_h, R_h. |
K_r |
Length- |
K_h |
Length- |
I_ini |
Length- |
R_ini |
Length- |
I_h_ini |
Length- |
R_h_ini |
Length- |
years |
Number of years to simulate (default 1). |
timestep |
"weekly" or "daily" (default "daily" – migration is typically a daily-cadence process). |
n_particles |
Number of particles for the stochastic run. |
n_threads |
Number of threads for parallel processing. |
obs_period |
Integer length (in days) of the observation window for
|
seasonal |
Optional length- |
... |
Additional scalar parameter overrides ( |
Per-patch parameters (K_r, K_h, I_ini, R_ini, I_h_ini, R_h_ini)
must be length-npop numeric vectors. All transmission and demographic
rates are scalars shared across patches (treat them as biological
constants, not place-specific). For per-patch beta_r[i] etc., edit
inst/odin/plague_stochastic_metapop.R directly.
plague_results tibble with population column running 1..npop.
Runs the stochastic rats + humans plague model
(inst/odin/plague_stochastic_humans.R) via dust2. The spatial rats-only
model was retired in April 2026 along with the dust v1 backend (see
archive/plague_stochastic.R). Re-adding spatial support means writing
a spatial+humans odin2 model first; this function will then grow a
npop argument again.
run_plague_model( scenario = "defaults", years = 10, timestep = c("weekly", "daily"), n_particles = 100, n_threads = 1, K_r = 2500, K_h = 5000, I_ini = 1, obs_period = 1L, ... )run_plague_model( scenario = "defaults", years = 10, timestep = c("weekly", "daily"), n_particles = 100, n_threads = 1, K_r = 2500, K_h = 5000, I_ini = 1, obs_period = 1L, ... )
scenario |
Parameter set name, file path, or list of parameters. |
years |
Number of years to simulate (default 10). |
timestep |
Time step resolution: "weekly" or "daily" (default "weekly"). |
n_particles |
Number of particles for the stochastic run. |
n_threads |
Number of threads for parallel processing. |
K_r |
Rat carrying capacity (used only when scenario does not set it). |
K_h |
Human carrying capacity (used only when scenario does not set it). |
I_ini |
Initial infected rats. |
obs_period |
Integer length (in days) of the observation window over
which |
... |
Additional parameters to override on top of the scenario. |
plague_results object.
Status bar module — server.
status_bar_server(id, lab_session)status_bar_server(id, lab_session)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
Invisibly, the moduleServer result.
Status bar module — UI.
status_bar_ui(id)status_bar_ui(id)
id |
Module namespace id. |
A shiny::div() styled as a sticky bottom bar.
Summarize outbreak metrics across replicates
summarize_outbreak_metrics(outbreak_metrics)summarize_outbreak_metrics(outbreak_metrics)
outbreak_metrics |
Output from calculate_outbreak_metrics |
Summary statistics
Summary method for plague_results
## S3 method for class 'plague_results' summary(object, ...)## S3 method for class 'plague_results' summary(object, ...)
object |
plague_results object |
... |
Additional arguments (ignored) |
Replaces each group's seasonal vector with group_seasonal[[g]] ^ alpha.
Lets the data inform the strength of seasonality — alpha = 0 collapses
seasonal to flat 1, alpha = 1 is the raw climate signal, intermediate
values dampen, larger values amplify. Removes alpha from unpacked pars.
with_alpha_seasonal(packer, group_seasonal)with_alpha_seasonal(packer, group_seasonal)
packer |
A grouped packer or wrapper around one. Must have |
group_seasonal |
Named list keyed by group; each entry is a numeric vector of per-day seasonal multipliers (length = number of simulation days for that outbreak). |
A packer with wrapped $unpack.
For each group, transforms the temperature time series into a per-day
multiplier on delta_R (carcass decay) using the Brière function
parameterised by fitted T_min, T_max, q_briere. The multiplier is
normalised so delta_R equals its scenario value at T_ref (default
17.8°C, the typical mid-range plague-permissive temperature) — i.e.
calibrating outside the Brière transformation reaches the unmodulated
scenario delta_R at the reference temperature.
with_briere_seasonal(packer, group_temp, T_ref = 17.8, cap = 1000)with_briere_seasonal(packer, group_temp, T_ref = 17.8, cap = 1000)
packer |
A grouped packer or wrapper around one. Must have |
group_temp |
Named list keyed by group; each entry is a numeric per-day temperature vector (°C, length = number of simulation days). |
T_ref |
Reference temperature at which the Brière multiplier equals
1 (i.e. the scenario |
cap |
Maximum allowed multiplier — clamps the transmission-off signal when params are out of range. Default 1000. |
If T_ref falls outside [T_min, T_max] (sampler exploring biologically
implausible cliffs), returns cap for all temperatures, which collapses
transmission and signals the proposal is bad. Per-temperature points
outside [T_min, T_max] are also capped. Removes T_min, T_max,
q_briere from unpacked pars.
A packer with wrapped $unpack.
For values that are fixed but vary per group (e.g. each outbreak's K_h,
K_r pinned to its Krauer-listed population), monty::monty_packer_grouped()'s
fixed argument doesn't help — it's global. This wrapper merges
group_fixed[[g]] into each group's pars during unpack.
with_per_group_fixed(packer, group_fixed)with_per_group_fixed(packer, group_fixed)
packer |
A grouped packer (from |
group_fixed |
Named list keyed by group; each entry is a named list of
parameters to splice into that group's unpacked pars. Groups not in
|
A packer with wrapped $unpack. Class preserved so dust2's grouped
filter accepts it.
Sampling R0 directly (rather than beta_r) gives a more interpretable prior
and a flatter posterior geometry. The conversion uses the disease-free
equilibrium derivation in the carcass model (see CLAUDE.md "R0 formula"):
beta_r = R0 * delta_R / ((1 - g_r) * (1 - exp(-rho))).
with_R0_to_beta_r(packer)with_R0_to_beta_r(packer)
packer |
A grouped packer or wrapper around one. Must have |
Requires g_r, rho, delta_R to already be present in each group's pars
(typically via the packer's global fixed argument). Removes R0 from the
unpacked pars after conversion.
A packer with wrapped $unpack.