| 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-10-22 11:33:51 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
Most of Krauer's catalogue is not fittable: 78 of the 130 records lack a
population (which lab_fit_assemble() pins K_h and K_r to), a
whole-day reporting window (which the D_h accumulator resets on), or
both. The explorer greys those out rather than offering them, so reaching
this error means an id arrived from a saved session or a script. Failing
here, naming the outbreak and the reason, beats an NA propagating into
dust2 and surfacing as an inscrutable sampler error.
cohort_check_fittable(cohort_ids, data = NULL)cohort_check_fittable(cohort_ids, data = NULL)
cohort_ids |
Character vector of |
data |
Long outbreaks tibble (defaults to bundled outbreaks). |
cohort_ids, invisibly, when every outbreak is fittable.
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.
Excludes parameters that aren't user-fittable in practice:
configurable_param_names()configurable_param_names()
tau (timestep — system, not biological)
seasonal_beta (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.
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).
Tile selection writes through to lab_session$cohort_ids, and the grid is
rebuilt from that vector, so a cohort changed anywhere else (loading a
saved fit, removing a chip in the Lab sidebar) is reflected here without a
second source of truth.
explore_server(id, lab_session, data = NULL)explore_server(id, lab_session, data = NULL)
id |
Module namespace id. |
lab_session |
A LabSession instance. |
data |
Long outbreaks tibble. Defaults to the bundled outbreaks_all. |
Invisibly, the moduleServer result.
Outbreak explorer — UI.
explore_ui(id)explore_ui(id)
id |
Module namespace id. |
A bslib::layout_sidebar() holding the filter rail, the card grid
and the detail card.
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.
Returns a shiny::shinyApp() object — call directly to launch
(lab_app()), or wrap with shiny::runApp() for explicit control over
host/port.
lab_app(saved_session = NULL, library_dir = NULL, restore_last = TRUE)lab_app(saved_session = NULL, library_dir = NULL, restore_last = TRUE)
saved_session |
Optional LabSession instance. If |
library_dir |
Directory backing the saved-fit library. |
restore_last |
Reopen the most recently saved fit (autosaves
included) when the app starts and no |
Requires the optional packages shiny, bslib, R6, ggplot2,
posterior, and DT. They live in Suggests; install with
install.packages(c("shiny", "bslib", "R6", "ggplot2", "posterior", "DT")).
A shiny::shinyApp() object.
Marginal posterior draws for several saved fits.
lab_compare_draws(slugs, library_dir = NULL, stages = NULL)lab_compare_draws(slugs, library_dir = NULL, stages = NULL)
slugs |
Character vector of library slugs (see |
library_dir |
Library directory. Defaults to |
stages |
Optional parallel vector of |
Long data frame with fit, variable, value. Zero rows when
nothing could be read.
A stochastic run stores pilot_samples alongside the production draws;
nothing else in the app can show them. This turns each such entry into two
comparable series so a run can be checked against its own pilot.
lab_compare_expand_stages(slugs, library_dir = NULL)lab_compare_expand_stages(slugs, library_dir = NULL)
slugs |
Character vector of library slugs. |
library_dir |
Library directory. Defaults to |
Entries with no separate pilot (a deterministic fit, where the pilot is the canonical posterior) contribute one series, never a duplicate.
List with slug and stage, suitable for lab_compare_draws().
Four. The categorical palette has four validated slots and hues are never cycled, so a fifth series would have no distinguishable colour.
lab_compare_max()lab_compare_max()
Integer.
Rebuilds each entry's dust2 setup (not serialised) and forward-simulates.
Entries whose setup cannot be rebuilt are skipped; their labels come back in
the skipped attribute so the caller can say which and why.
lab_compare_trajectories( slugs, library_dir = NULL, n_draws = 40L, stages = NULL )lab_compare_trajectories( slugs, library_dir = NULL, n_draws = 40L, stages = NULL )
slugs |
Character vector of library slugs. |
library_dir |
Library directory. Defaults to |
n_draws |
Posterior draws per fit. Default 40 — lower than the hero's, since several fans share one panel. |
stages |
Optional parallel vector of |
Long data frame with fit, draw, group, time, mu, carrying
a skipped attribute (character vector of labels).
Returns a list containing the posterior model, sampler, packer, fixed
parameters, and metadata. Pass to lab_fit_run() (or call
monty::monty_sample() directly with posterior + sampler).
lab_fit_assemble(lab_session, data = NULL)lab_fit_assemble(lab_session, data = NULL)
lab_session |
A LabSession whose |
data |
Long outbreaks tibble (defaults to bundled outbreaks). |
Branches on length(cohort_ids):
1: standard single-outbreak setup, plain monty_packer.
>1: grouped fit with monty::monty_packer_grouped() and
with_per_group_fixed() splicing per-outbreak populations into K_h
and K_r when those aren't in the fitted set.
Smoke-tests the unfilter at the prior medians so misconfigured fits
fail at assembly time, not deep inside monty_sample().
Named list: posterior, sampler, packer, unfilter,
fixed_pars, data, cohort_ids, fitted_names.
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, deterministic = TRUE, seed = NULL )lab_fit_forward_sim( setup, samples, n_draws = 100, t_grid = NULL, deterministic = TRUE, seed = 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
|
deterministic |
Run the model deterministically ( Deterministic runs show parameter uncertainty alone: each draw traces
the expected epidemic for its parameters. That matches the deterministic
unfilter the pilot was fitted with, so it is the honest picture of a
pilot. A stochastic fit was fitted against a particle filter, where
demographic stochasticity is part of the likelihood — extinction of a
small infected seed, the timing jitter of an epizootic — and a
deterministic fan hides exactly the variation that fit was paying for.
Pass |
seed |
Optional integer seed for the stochastic path, for
reproducible fans. Ignored when |
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
|
The second half of the two-stage workflow: a deterministic pilot tunes the proposal, then a short stochastic production run starts from it.
A monty_samples object with the warmup iterations removed.
Called by the app after every successful fit. Generates a name from the cohort, scenario, mode and timestamp, and writes a new versioned entry — autosaves never replace one another, so re-running the same configuration leaves both runs in the library to compare.
lab_library_autosave(lab_session, library_dir = NULL)lab_library_autosave(lab_session, library_dir = NULL)
lab_session |
A LabSession instance with a completed fit. |
library_dir |
Directory to write to. Defaults to |
Failures are swallowed deliberately: a full disk or an unwritable library directory should never turn a successful fit into an error the user sees. The condition message is returned invisibly for logging.
Invisibly, the metadata list, or NULL if the save failed.
Delete a saved fit from the library.
lab_library_delete(slug, library_dir = NULL)lab_library_delete(slug, library_dir = NULL)
slug |
Slug of the saved fit. |
library_dir |
Directory to delete from. Defaults to |
Invisibly, TRUE if anything was removed.
Defaults to a per-user data directory so saved fits survive across R
sessions and package reinstalls. Override with the yersinia.library_dir
option or by passing library_dir to lab_app().
lab_library_dir(create = TRUE)lab_library_dir(create = TRUE)
create |
Create the directory if it doesn't exist. Default |
Path to the library directory.
Strictly most-recent by timestamp, autosaves included — unlike
lab_library_list()'s display order, which floats named entries to the
top. Used to restore the last fit when the app opens.
lab_library_last(library_dir = NULL)lab_library_last(library_dir = NULL)
library_dir |
Directory to scan. Defaults to |
A slug, or NULL when the library is empty.
Reads only the .meta.rds sidecars, so listing stays fast no matter how
large the stored posteriors get. Entries whose sidecar is unreadable are
skipped rather than aborting the listing.
lab_library_list(library_dir = NULL)lab_library_list(library_dir = NULL)
library_dir |
Directory to scan. Defaults to |
A data frame, one row per saved fit, newest first. Zero rows (with the full set of columns) when the library is empty.
Restores cohort, model config, priors, and the stored fit, then rebuilds
the dust2 setup that was deliberately not serialised. If the rebuild fails
(e.g. the bundled outbreak data changed under the saved cohort), the rest
of the state is still applied and setup is left NULL — the fit is
viewable, it just can't be promoted to production without a refit.
lab_library_load(lab_session, slug, library_dir = NULL)lab_library_load(lab_session, slug, library_dir = NULL)
lab_session |
A LabSession to mutate in place. |
slug |
Slug of the saved fit (from |
library_dir |
Directory to read from. Defaults to |
Invisibly, lab_session.
Clears the autosave flag and renames the entry, so a run worth keeping stops looking like scratch and sorts with your named fits. The payload is re-read and re-written under the new slug; the old pair is removed.
lab_library_promote(slug, name, library_dir = NULL)lab_library_promote(slug, name, library_dir = NULL)
slug |
Slug of the autosaved entry. |
name |
New human-readable name. |
library_dir |
Directory to operate in. Defaults to |
Invisibly, the new metadata list.
Writes lab_session$to_list() plus a metadata sidecar.
lab_library_save(lab_session, name, library_dir = NULL, autosaved = FALSE)lab_library_save(lab_session, name, library_dir = NULL, autosaved = FALSE)
lab_session |
A LabSession instance. |
name |
Human-readable name for this fit. Used to derive the filename. |
library_dir |
Directory to write to. Defaults to |
autosaved |
Mark this entry as automatically saved. Autosaves are flagged in the metadata so the app can filter them out of the picker, and are never overwritten. |
Two collision policies, because manual and automatic saves want opposite
things. A manual save overwrites any existing entry with the same slug,
so iterating on one named fit updates it in place instead of accumulating
near-duplicates. An autosave (autosaved = TRUE) never overwrites: its
generated name already carries a timestamp, and a numeric suffix is added
in the unlikely event two runs land in the same second. Every run is
therefore recoverable.
Invisibly, the metadata list that was written.
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
Maps the nine hand-assigned string ids the app used before outbreaks were
rekeyed on Krauer's catalogue id to their current ids, so saved sessions
and fit-library entries keep resolving. Consumed by
outbreak_resolve_id().
outbreak_aliasesoutbreak_aliases
A tibble with 9 rows and 2 variables:
The old string, e.g. "Barcelona_1490".
The current id.
"Barcelona_1490" maps to the outbreak now labelled Barcelona
1489. The year really does change: Krauer's record begins 1489-11-05, and
the old hand-transcribed file dropped its first 125 days. The other eight
legacy ids matched their natural spelling exactly.
Display labels for outbreak identifiers.
outbreak_label(ids, data = NULL)outbreak_label(ids, data = NULL)
ids |
Character vector of outbreak identifiers, current or legacy. |
data |
Long-format outbreaks tibble. Defaults to outbreaks_all. |
Character vector of labels; unresolvable ids are returned as-is.
outbreak_label(c("18", "69"))outbreak_label(c("18", "69"))
Current ids pass through untouched. Legacy string ids (the nine the app used before outbreaks were rekeyed on Krauer's catalogue id) are translated via outbreak_aliases. Ids matching neither are returned unchanged, so the caller's own "no data for this outbreak" error is what the user sees rather than one from here.
outbreak_resolve_id(ids, data = NULL)outbreak_resolve_id(ids, data = NULL)
ids |
Character vector of outbreak identifiers, current or legacy. |
data |
Long-format outbreaks tibble used to decide what counts as a current id. Defaults to the bundled outbreaks_all. |
Note that "Barcelona_1490" resolves to the outbreak now labelled
Barcelona 1489. That is not an off-by-one: Krauer's record begins
1489-11-05, and the old hand-transcribed file dropped its first 125 days,
which is what made it look like a 1490 outbreak.
Character vector the same length as ids.
outbreak_resolve_id(c("London_1563", "31"))outbreak_resolve_id(c("London_1563", "31"))
Aggregates the long-format outbreaks_all dataset to one row per outbreak, with the metadata, summary statistics and screening flags the explorer filters, sorts, maps and flags on.
outbreak_summary(data = NULL)outbreak_summary(data = NULL)
data |
Long-format outbreaks tibble. Defaults to the package's bundled outbreaks_all dataset. |
The returned columns beyond the obvious ones:
labelDisplay string ("Barcelona 1489"), disambiguated by
start month where place and year collide. For reading, never joining —
join on outbreak_id.
fittableWhether the app can fit this record at all. Requires
both a population (lab_fit_assemble() pins K_h and K_r to it) and
a fixed reporting window (the D_h accumulator resets every
obs_period steps, which monthly records have no whole-day value for).
52 of the 130 qualify.
unfit_reasonWhich gate failed, in a phrase fit to show the
user. NA when fittable.
attack_rateRecorded deaths divided by population.
attack_flag"impossible" above 100% (Klaipeda 1710 records
9,797 deaths in a population of 4,000), "high" above 35% — the
threshold at which a closed-population model stops being able to
reproduce the series. It is a warning, not a filter: Prague 1713 and
Eyam 1665 both sit at 37% and are in the curated nine.
peak_doyDay of year (1-366) on which the largest observed count falls – the outbreak's seasonal phase. Correlates with latitude at r = 0.64 across the catalogue: north of 52 N the mean peak is day 242, south of 45 N it is day 180.
type"plague mortality" or "all-cause mortality". 54 of the
130 are all-cause burials, for which lambda_baseline is carrying a
different quantity — see the fitting notes in CLAUDE.md.
A tibble with one row per outbreak_id, sorted by year.
outbreak_summary() # Just the ones that can actually be fitted: subset(outbreak_summary(), fittable)outbreak_summary() # Just the ones that can actually be fitted: subset(outbreak_summary(), fittable)
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) )
All 130 outbreaks from Fabienne Krauer's plague-season v3.2 dataset (84 places, 22 countries, 1348–1878), in the same long format and with the same key scheme as outbreaks, which is a nine-record subset of it. This is the dataset the outbreak explorer works from.
outbreaks_alloutbreaks_all
A tibble in long format, one row per observation:
Krauer's catalogue id as a character string. Unique
by construction. Not place_startyear, which is not:
Krauer 5 and 6 are both Alexandria 1840, and 7 and 8 are both
Alexandria 1842.
Display string, e.g. "Barcelona 1489", disambiguated
by start month where place and year collide. For reading, not joining.
Place and time metadata.
Total population, or NA. Missing for 61 of the
130 records.
Day index from the outbreak start. For daily series the day itself; for weekly and biweekly the day the reporting window closes; for monthly the day it opens.
Count for the window ending at day. NA means
no observation and must never be read as an observed zero.
Reporting-window length in days (1, 7 or 14), or
NA for monthly series, which have no whole-day window.
Provenance, and the fields the seasonal-forcing work joins on.
"plague mortality" or "all-cause mortality";
54 of the 130 are all-cause burials.
Krauer's own completeness flag, and whether the series was transcribed from a table or digitised from a graph.
Whether the app can fit this record: needs both a population and a whole-day reporting window. 52 of the 130 qualify.
Which requirement failed; NA when fittable.
Recorded deaths, and that over population.
"high" above 35% and "impossible"
above 100%. A warning, not a screen — Prague 1713 and Eyam 1665 sit
at 37% and are both in the curated outbreaks subset.
Krauer, F. plague-season v3.2 (input/rawdata.csv), vendored
at data-raw/krauer-plague-season-v3_2.csv.
outbreak_summary() for the one-row-per-outbreak view,
outbreak_resolve_id() for translating pre-2026-09 identifiers.
# How much of the catalogue is actually modellable, and why not: table(outbreak_summary()$unfit_reason, useNA = "ifany")# How much of the catalogue is actually modellable, and why not: table(outbreak_summary()$unfit_reason, useNA = "ifany")
Renders beta_h as "Carcass-to-human transmission (beta_h)", with the
symbol in Greek where the literature uses Greek. Grouped parameter names
from monty::monty_packer_grouped() (param<group>) keep their group,
appended after a separator, so per-outbreak panels stay distinguishable.
param_label(x, with_symbol = TRUE)param_label(x, with_symbol = TRUE)
x |
Character vector of parameter names. |
with_symbol |
Include the symbol in parentheses. Default |
Unknown names pass through unchanged – a parameter added to the model still plots, it just isn't prettified yet.
Character vector of labels, same length as x.
param_label(c("beta_h", "kappa", "beta_h<Eyam_1665>"))param_label(c("beta_h", "kappa", "beta_h<Eyam_1665>"))
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_beta = 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_beta = 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_beta |
Optional per-step thermal multiplier on |
... |
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, library_dir = NULL)status_bar_server(id, lab_session, library_dir = NULL)
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) |
Returns a per-temperature multiplier on transmission, equal to 1 at
T_opt and falling to zero at T_opt - hw_cold and T_opt + hw_hot.
This is simplified Brière-2 reparameterised so the peak and the distances
to the two zeros are the parameters:
thermal_response(temperature, T_opt, hw_cold, hw_hot, floor = 1e-08)thermal_response(temperature, T_opt, hw_cold, hw_hot, floor = 1e-08)
temperature |
Numeric vector of temperatures (degrees Celsius). |
T_opt |
Temperature at which transmission peaks. |
hw_cold |
Degrees from |
hw_hot |
Degrees from |
floor |
Smallest value returned, applied outside the zeros and for degenerate parameters, so a bad MCMC proposal lowers the likelihood instead of producing zero transmission or an error. Default 1e-8. |
Equivalent textbook coordinates: T_min = T_opt - hw_cold,
T_max = T_opt + hw_hot, exponent q = hw_hot / hw_cold.
Numeric vector in [floor, 1], one per input temperature.
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.
For each group, turns that group's daily temperature series into a per-day
multiplier on beta_r and beta_h via the odin model's seasonal_beta
input, using thermal_response() with fitted T_opt, hw_cold, hw_hot.
Carcass decay is not forced.
with_thermal_beta(packer, group_temp)with_thermal_beta(packer, group_temp)
packer |
A grouped packer or wrapper around one. Must carry |
group_temp |
Named list keyed by group; each entry is a numeric per-day temperature vector (degrees C, length = number of simulation days). |
Bad proposals (non-positive widths) yield a floored response rather than an error, which drives the likelihood down and pushes the sampler away.
A packer with wrapped $unpack.