Clean up workflow ergonomics, defer flat_prior-only inits, drop sero for now#52
Merged
Conversation
…series The date-alignment check in prepare_observation_data() assumed any already- classed epiwave_timeseries object stores its dates in a flat $date column, but epiwave_greta_timeseries objects (e.g. the IHR-from-CHR pattern used for hospitalisation proportions in tests/test_workflow/testing.R and test2.R, built via create_epiwave_greta_timeseries()) are a list wrapping a greta array alongside the date tibble, with dates nested at $timeseries$date instead. The check was comparing as.Date(NULL) against target_infection_dates and always failing with "`proportion_infections` dates must match `target_infection_dates`". Found by actually running tests/test_workflow/testing.R's basic (cases + hospitalisations) fit against real data in ../data -- this is a core supported pattern (proportion driven by a greta array), not an edge case, and would have broken on the first real multi-stream fit using it. Verified: the same fit now runs to completion against real data.
Adds as_epiwave_timeseries(data): a new exported helper that takes a plain
data.frame/tibble with date/value columns (partial date coverage is fine)
and classes it as epiwave_fixed_timeseries -- replacing the fragile,
unvalidated class(x) <- c('epiwave_fixed_timeseries', 'epiwave_timeseries',
class(x)) pattern users previously had to hand-write for every observation
stream. Already-classed epiwave_timeseries objects and plain numeric
values/vectors pass through unchanged, so it's safe to call unconditionally.
prepare_observation_data() now calls this automatically on timeseries_data
and size_vec (mirroring the auto-coercion that delay_from_infection/
proportion_infections already had), so users can pass a plain data.frame
straight into define_observation_data()/define_sero_data() without ever
manually setting a class themselves. This directly targets the messy
"before define_observation_data()" data prep in tests/test_workflow/
testing.R, where every observation stream needed its own hand-rolled
class(x) <- ... line.
Verified: a raw unclassed data.frame produces identical output to the
equivalent pre-classed object; a data.frame missing the required columns
errors with a clear message instead of failing deep inside as_matrix().
define_observation_data() now accepts a plain date/value data.frame directly
(via as_epiwave_timeseries(), added earlier on this branch), so the
hand-rolled class(notif_dat) <- c('epiwave_fixed_timeseries',
'epiwave_timeseries', class(notif_dat)) lines are no longer needed.
Verified: re-ran testing.R's basic (cases + hospitalisations) fit against
real data in ../data with the classing removed -- notif_dat/hosp_dat stay
plain tbl_df objects throughout, and the fit runs to completion identically.
testing.R/test2.R depend on local, not-synced data (../data, simdata/) that isn't in the repo -- useful as personal real-data scripts, but not runnable by anyone else who clones it, and don't exercise the multi-jurisdiction path at all (the package has never had a multi-jurisdiction reprex before this branch). Adds two fully self-contained scripts (fabricated data, no external dependencies) that demonstrate the current API end to end: - single_jurisdiction_workflow.R: raw data.frames passed straight to define_observation_data() (no manual class(x) <- c(...)), delay distributions built via epiwave.params (discrete_pmf, combined with `+`), a greta-array-derived proportion (IHR-from-CHR), and an "advanced" section showing a time-varying discrete_pmf_series delay. - multi_jurisdiction_workflow.R: two jurisdictions with deliberately staggered/non-identical date coverage, dow_model = TRUE to exercise hierarchical DOW pooling, and an explicit alignment check confirming a date only one jurisdiction has data for stays correctly non-NA/NA in the right columns. Both run to completion (small MCMC settings for speed) with their stopifnot() checks passing.
Captures why the seroprevalence pathway (define_sero_data(), create_small_sero_model(), the discrete_weights widening in new_convolution_matrix()/evaluate()) is being pulled back out of this MVP, the recommended architecture for re-integrating it (merge into define_observation_data()/create_observation_model() rather than parallel functions with if-based dispatch), the real bugs found in the original sero code during this work, and what validation is still needed before it ships. dev/ is excluded from R CMD build via .Rbuildignore, matching the existing pattern for .claude/.positai.
Per discussion: don't ship if-statement dispatch for a not-yet-validated feature. Removes define_sero_data(), create_small_sero_model(), the total_pop/size_vec plumbing through prepare_observation_data() and stack_jurisdictions()'s stack_stream(), the model_fn dispatch in fit_waves(), and the discrete_weights/discrete_weights_series widening in new_convolution_matrix()/evaluate() (which had no other consumer once sero is removed). inits_by_jurisdiction()'s mean_step() helper simplifies back to a plain mean(x) call accordingly. The design for how to re-integrate this cleanly (merge into define_observation_data()/create_observation_model() rather than parallel functions with if-based dispatch, as agreed before removing it), the real bugs found in the original sero code, and what's still needed before it ships are captured in dev/sero-integration-notes.md (previous commit). fit_waves()/stack_jurisdictions()/prepare_observation_data() now only handle the validated cases/hospitalisations pathway.
testing.R's sero fit block called define_sero_data(), which no longer exists after the previous commit -- it was already non-functional (sero_dat/ sero_size_mat/sero_conversion were never assigned), so this just removes dead code rather than breaking anything that worked. Left short pointer comments to dev/sero-integration-notes.md at both spots (data prep, and where a third define_observation_data() block would go) for whoever picks this back up. devtools::document() drops create_small_sero_model/define_sero_data exports and man pages, and the now-unused importFrom(greta, binomial). Verified after all sero-removal changes: devtools::check() is 0 errors, 0 warnings (same 3 pre-existing NOTEs); testing.R's basic fit still runs against real data in ../data; both single_jurisdiction_workflow.R and multi_jurisdiction_workflow.R still run to completion with their stopifnot() checks passing.
Two changes, implemented together since they touch the same functions: 1. stack_jurisdictions() becomes the explicit, user-facing combining step, taking named ... args (mirroring define_observation_model()'s own cases=/hospitalisations= pattern) instead of a pre-built named list -- e.g. stack_jurisdictions(VIC = obs_vic, NSW = obs_nsw). A single jurisdiction needs no combining step at all: define_observation_model()'s output (now stamped with class epiwave_observation_model) goes straight to fit_waves(), which detects whether it already received a stacked object (class epiwave_stacked_observations) or a raw single-jurisdiction one and wraps the latter internally (auto-labelled, since there's no jurisdiction identity to preserve for n=1). This removes the setNames(list(...), jurisdictions) step that was previously required even for a single jurisdiction. 2. GAM-based initial values (inits_by_jurisdiction(), and the cross-stream union/mean that combines them) are only ever consumed by infection_model_type = 'flat_prior' -- every GP-based infection model ignores them entirely (confirmed in create_infection_timeseries(): observable_infection is referenced only inside the flat_prior branch). Previously this ran unconditionally for every stream x jurisdiction during data prep, running an mgcv::gam() fit that was silently discarded whenever a GP model was used instead. prepare_observation_data() no longer computes inits at all (keeps delays in its output instead, needed later); define_observation_model()/stack_jurisdictions() no longer eagerly combine them either. A new compute_flat_prior_inits() does this work on the fully-stacked object, called lazily by fit_waves() only inside its flat_prior branch, before create_infection_timeseries() (which needs the result) and reused for the MCMC initial values. Verified: a standalone check exercising both a single-jurisdiction fit (flat_prior AND gp_growth_rate, confirming the GP path never triggers inits computation) and a two-jurisdiction stack_jurisdictions() combination (plus the unnamed-args error path) all pass. devtools::check() is 0 errors, 0 warnings. testing.R's real-data fit (../data) verified against both flat_prior and gp_growth_rate. Both synthetic workflow scripts updated to the new API and re-verified end to end.
compute_flat_prior_inits() was reading prop_mat from the stacked object, but stack_stream() applies DOW correction to prop_mat before returning it -- so for any stream with dow_model = TRUE, inits_by_jurisdiction() was being handed a greta array (the DOW-corrected proportion) instead of a plain numeric one. It didn't error (there's already a greta_array branch in inits_by_jurisdiction()), but that branch draws from the *prior* predictive distribution of the DOW effect via greta::calculate(nsim = 100) -- no MCMC has run yet at this point -- making the "deterministic smoothed guess" inits computation silently stochastic, and adding an unnecessary 100-draw simulation per DOW-modelled stream x jurisdiction whenever flat_prior is used. Checked against the original, pre-refactor code to confirm this wasn't just a style difference: it explicitly computed inits from prop_mat *before* applying DOW correction, deliberately. My first refactor (PR #51) preserved this ordering (DOW correction lived in stack_jurisdictions(), which ran after prepare_observation_data()'s inits computation); deferring inits to run after stacking (previous commit) inverted it by accident. Fix: stack_stream() now also returns prop_mat_raw (captured before DOW correction is applied), and compute_flat_prior_inits() reads that instead of prop_mat. create_observation_model() is unaffected -- it still reads the DOW-corrected prop_mat for the actual likelihood, as it should. Verified: prop_mat_raw is a plain numeric column (not a greta_array) even when dow_model = TRUE, so inits_by_jurisdiction() no longer needs greta::calculate() at all for this case. Re-ran both synthetic workflow scripts and testing.R's real-data fit (flat_prior + dow_model = TRUE, the exact previously-affected combination) end to end. devtools::check() is still 0 errors, 0 warnings.
It's called from exactly one place (define_observation_model()) and never appears in any workflow script -- it's implementation detail, the same category as stack_stream()/compute_flat_prior_inits(), which are already @nord. Exporting it overstated it as public API and was contributing to a "why are there three equally-important functions here" feeling, when really only define_observation_data()/define_observation_model() (plus stack_jurisdictions()/fit_waves()) are meant to be called directly. No logic changed -- devtools::check() still 0 errors/0 warnings, both synthetic workflow scripts re-verified end to end.
Same situation as prepare_observation_data(): called from exactly one place (compute_flat_prior_inits()), never from a workflow script. Pure internal machinery. No logic changed -- devtools::check() still 0 errors/0 warnings, both synthetic workflow scripts re-verified end to end.
… simplify prop coercion create_epiwave_fixed_timeseries() had zero call sites anywhere in the package or workflow scripts, and did nothing as_epiwave_timeseries() doesn't already do better (it takes a real date/value table directly, rather than requiring separate dates + value arguments). create_epiwave_timeseries() had exactly one remaining caller, in prepare_observation_data()'s proportion_infections coercion -- but that coercion step was itself redundant: as_matrix() already has an as_matrix.numeric method that recycles a scalar or validates a same-length vector against target_infection_dates, identically to what going through create_epiwave_timeseries() first produced (verified: identical() output both ways, for both the scalar and vector cases). prepare_observation_data() now only branches on already-classed epiwave_timeseries objects (for date validation); a bare numeric prop passes straight to as_matrix(), which dispatches correctly on its own. create_epiwave_timeseries.R now holds only create_epiwave_greta_timeseries() (kept in epiwave rather than epiwave.params specifically because it depends on greta, which epiwave.params should not) and as_epiwave_timeseries(). Verified: both synthetic workflow scripts and testing.R's real-data fit produce identical output to before these changes (e.g. single-jurisdiction posterior median range unchanged at 0-7185). devtools::check() still 0 errors, 0 warnings.
…-export as_epiwave_timeseries() Renamed to pair naturally with as_epiwave_timeseries() now that the two create_*/dead functions removed from this file are gone -- as_* matches the coercion-style naming already used elsewhere (epiwave.params::as_discrete_pmf(), as_matrix()). The class it produces is renamed to match: epiwave_greta_timeseries -> greta_timeseries (still inherits from epiwave_timeseries, so existing inherits(x, 'epiwave_timeseries') checks that need to catch both plain and greta-backed timeseries objects continue to work unchanged). as_matrix.epiwave_greta_timeseries() is renamed to as_matrix.greta_timeseries() to match, and prepare_observation_data()'s inherits() check on proportion_infections updated accordingly. as_epiwave_timeseries() becomes internal (@nord): checked and confirmed no workflow script actually calls it by name -- users get its behaviour indirectly (passing a raw data.frame to define_observation_data(), which coerces internally). It remains available as an unexported helper rather than being deleted, since it's still the one general-purpose coercion path, just not meant as public API for now. Verified: class(ihr) is now `greta_timeseries epiwave_timeseries list`, dispatch to as_matrix.greta_timeseries() works correctly, both synthetic workflow scripts produce identical output to before the rename (e.g. single-jurisdiction posterior median range unchanged at 0-7185), and testing.R's real-data fit runs to completion. devtools::check() still 0 errors, 0 warnings.
Two minimal examples -- single jurisdiction, and multiple jurisdictions combined via stack_jurisdictions() -- mirroring tests/test_workflow/single_jurisdiction_workflow.R and multi_jurisdiction_workflow.R, trimmed to the essential shape for a README. Both verified to run to completion exactly as written before adding them; chunks are eval = FALSE in README.Rmd since running them needs a full greta/Python setup, which shouldn't be a precondition for a routine README re-knit. README.Rmd had been missing from the repo (only README.md existed) despite a pre-commit hook expecting both to move together; it's restored here from an old, partially-filled scaffold (its Installation section still suggested pak::pak() and its Example section was a never-filled-in placeholder). Updated Installation to match the remotes::install_github() already decided in README.md, and replaced the Example placeholder with the new Usage section. Left the empty Citation section and placeholder Support text untouched rather than inventing content, and left the informal draft paragraph after the srr-tags chunk untouched (out of scope for this change). Also gitignores README.html, a preview-render artifact that isn't source.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Follow-up to #51, addressing review feedback on that branch before it merges further, plus general polish for this week's MVP.
Ergonomics — no more manual jurisdiction naming/listing.
stack_jurisdictions()is now the explicit, user-facing combining step, taking named...args (mirroringdefine_observation_model()'s owncases=/hospitalisations=pattern) instead of a pre-builtsetNames(list(...), jurisdictions). A single jurisdiction needs no combining step at all —define_observation_model()'s output goes straight tofit_waves(), which detects whether it received an already-stacked object or a raw single-jurisdiction one (via class) and wraps the latter internally.Deferred, flat_prior-only inits. GAM-based initial values (
inits_by_jurisdiction()) were running unconditionally for every stream × jurisdiction during data prep, even though GP-based infection models never reference them at all (confirmed increate_infection_timeseries()—observable_infectionis only used inside theflat_priorbranch). A newcompute_flat_prior_inits()computes this lazily, only insidefit_waves()'sflat_priorbranch, on the fully stacked object.A real bug found and fixed along the way:
compute_flat_prior_inits()initially read the DOW-correctedprop_matfrom the stacked object, but the original (pre-refactor) code deliberately computed inits from the proportion before DOW correction — using the corrected version made the "deterministic smoothed guess" depend on a prior-predictive draw of the DOW effect instead. Fixed by havingstack_stream()retainprop_mat_raw(pre-correction) specifically for this.Seroprevalence pathway removed from this MVP (not merged, not abandoned):
define_sero_data(),create_small_sero_model(), and thediscrete_weightswidening innew_convolution_matrix()/evaluate()are pulled back out — the pathway was only ever validated against fabricated data, and we didn't wantif-based dispatch for a not-yet-validated feature sitting in the MVP's core functions. Design notes for reintroducing it cleanly (merge intodefine_observation_data()/create_observation_model()rather than parallel functions) are indev/sero-integration-notes.md(excluded from the package build).API surface cleanup.
prepare_observation_data()andinits_by_jurisdiction()are now internal (@noRd) — both were only ever called from one place, never from a workflow script.create_epiwave_timeseries()/create_epiwave_fixed_timeseries()are deleted entirely: the former had exactly one caller whose jobas_matrix()'s existingnumericdispatch already did on its own (verified identical output), the latter had zero callers anywhere.create_epiwave_greta_timeseries()is renamed toas_greta_timeseries()(class renamedepiwave_greta_timeseries→greta_timeseriesto match) to pair naturally with the newas_epiwave_timeseries()auto-coercion helper, which replaces the need for users to hand-writeclass(x) <- c(...)on their raw observation data.New self-contained workflow scripts and README usage examples.
tests/test_workflow/{single,multi}_jurisdiction_workflow.Rare fully self-contained (fabricated data, no external files), unliketesting.R/test2.Rwhich depend on local not-synced data. The multi-jurisdiction one deliberately uses staggered/non-identical date coverage between jurisdictions to exercise the master-date-axis alignment.README.Rmd/README.mdgained a Usage section with trimmed versions of both.Test plan
devtools::check()— 0 errors, 0 warnings throughout (same 3 pre-existing NOTEs unrelated to this branch)dow_model = TRUE+flat_priorcombination that exposed the DOW/inits bugtesting.R's real-data fit (../data) re-verified against bothflat_priorandgp_growth_rateafter every changestack_jurisdictions()'s named-args API, the unnamed-args error path,fit_waves()correctly routing both a raw single-jurisdiction object and an already-stacked multi-jurisdiction object, and confirminggp_growth_ratenever triggers the GAM-based inits computation at allclass(ihr)confirmedgreta_timeseries epiwave_timeseries listpost-rename, withas_matrix.greta_timeseries()dispatch working🤖 Generated with Claude Code