Two station static model: Stan file, data, data prep, implementation by mdodrill-usgs · Pull Request #469 · DOI-USGS/streamMetabolizer

mdodrill-usgs · 2026-07-13T19:40:50Z

Makes progress on #441: Implements the two-station (VFTS) metabolism model as a new model type
in streamMetabolizer. All work targets the dev/two-station integration
branch.

Changes

inst/models/b2_np_oi_tr_plrckm.stan — static Stan model (fixed
syntax, lognormal K600, GPP/ER priors, package naming conventions)
R/mm_data.R — new two-station columns (DO.obs.up, DO.sat.up,
DO.obs.down, DO.sat.down, travel.time)
R/mm_validate_data.R — fix pre-existing travel.time grep bug in
na_times check
R/mm_ts_prep_data.R — data reshape with upstream lag shift
R/metab_2station.R — metab_2station() implementation +
predict_DO.metab_2station()
R/mm_name.R, R/mm_valid_names.R — register bayes_2station type
R/mm_parse_name.R — extend type map for b2_ prefix
R/specs.R — bayes_2station specs case; prior defaults match
one-station convention (K600: meanlog=2.484907, sdlog=1.0)
R/metab.R — dispatch to metab_2station()
R/data.R, data/two_station_example.rda — 30-day example dataset
(Bishop et al. 2026, lom3.70066)
DESCRIPTION — add stats, utils, parallel to Imports

Tests

472 passing tests including a real Stan fit against example data.
0 errors | 3 warnings (pre-existing) | 6 notes (5 pre-existing,
1 resolved in this PR)
Test timing: ~192 sec total, metab_2station adds ~20 sec

Known gaps (deferred, noted in code)

predict_DO(use_saved=FALSE) not implemented
plot_DO_preds.R and helper-rmse_DO.R hard-code one-station column
names

References

Bishop et al. 2026. A simplified two-station approach for modeling
metabolism in dam tailwaters subject to diel flow variation.
L&O Methods. https://doi.org/10.1002/lom3.70066

Port VFTS Stan model from LOM paper to inst/models/ - Fix deprecated postfix array syntax (10 declarations; Stan 2.33+ requirement) - Switch K600 prior from hard-coded normal to lognormal with data-block parameters - Rename variables to package conventions (n_obs, n_days, light, temp_water, travel_time) - Verified clean parse

- Add DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, travel.time to mm_data() registry - Create metab_2station() constructor stub with mm_validate_data() integration - Two-station-specific checks: travel.time positivity, units sanity (<1 day), lead-in data - Fix mm_validate_data() na_times grep incorrectly matching travel.time alongside solar.time - Add test-metab_2station.R (7 assertions); all pre-existing tests still pass

- Implement mm_ts_prep_data(data, specs=NULL) reshaping long-format data to Stan list - Upstream DO/DO.sat shifted per-row by lag = round(travel.time / timestep_days) - Lead-in rows (first max(lag) rows) excluded from output matrices - Pivots 8 variables to [n_obs x n_days] matrices following prepdata_bayes() pattern - K600 prior params passed through from specs (defaults: log(3.48) / 0.5) - Units stripped before return - Extend test-metab_2station.R to 15 blocks / 39 assertions; no regressions

- Extend mm_parse_name() type map: b2 = 'bayes_2station' (handles two-character prefix correctly) - Add predict_DO.metab_2station() stub returning DO.obs.down / DO.mod.down column names - Flag plot_DO_preds.R and helper-rmse_DO.R as needing update - Extend test-metab_2station.R to 12 assertions; no regressions in test-mm_name.R

- Rename GPP/ER/k600 to GPP_daily/ER_daily/K600_daily to match format_mcmc_mat_nosplit() lookup table - Add GPP_daily_mu/sigma and ER_daily_mu/sigma to data block - Add GPP_daily and ER_daily normal priors in model block - Group all three priors before likelihood statement

- Register bayes_2station in mm_name(), mm_valid_names() - Add bayes_2station case to specs(): params_in/out, split_dates=FALSE, engine='stan' - K600 prior defaults match one-station convention (meanlog=2.484907, sdlog=1.0) - GPP/ER prior defaults match one-station defaults (GPP: mu=3.1 sd=6.0, ER: mu=-7.1 sd=7.1) - Add stats, utils, parallel to DESCRIPTION Imports (pre-existing gap) - Add bayes_2station dispatch to metab() - Full metab_2station() implementation: validation, data prep, rstan::stan() call, MCMC formatting via format_mcmc_mat_nosplit(), metab_model S4 construction - Add metab_2station class - Implement predict_DO.metab_2station() returning DO.obs.down / DO.mod.down - Add two_station_example dataset (30 days, VFTS-2, package column conventions) - Add R/data.R with full citation: Bishop et al. 2026 and ScienceBase link - passing tests including real Stan fit; no regressions in mm_name suite - Known gaps: no Stan compile caching (deferred), predict_DO use_saved=FALSE not implemented, plot_DO_preds.R/helper-rmse_DO.R one-station column names deferred

- Add utils::globalVariables for ., metab_50pct, DO.mod.down - Add @importFrom utils modifyList - Replace require(rstan) with requireNamespace() pattern

mdodrill-usgs · 2026-07-13T20:11:18Z

Demo

Ran the two-station model against the bundled two_station_example dataset to confirm the full fit → predict pipeline works.

Dataset

?two_station_example
two_station_example

dim: 2904 x 9
names: solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, depth, temp.water, travel.time

Model fit

specs <- specs(mm_name('bayes_2station'),
               n_chains = 1, burnin_steps = 250, saved_steps = 250)
fit <- metab(specs, data = two_station_example)

Fast demo settings (1 chain, 250 warmup + 250 saved steps) to confirm the pipeline runs.

`predict_metab()` — daily estimates (first 5 days)

date	GPP	GPP.lower	GPP.upper	ER	ER.lower	ER.upper	K600	K600.lower	K600.upper
2011-07-31	16.88	15.28	18.23	-31.63	-38.54	-24.31	23.26	15.93	30.39
2011-08-01	18.31	16.68	20.06	-30.07	-37.53	-24.80	20.72	15.65	28.43
2011-08-02	19.74	18.30	21.15	-33.00	-39.26	-27.05	24.28	17.99	30.87
2011-08-03	18.62	17.31	19.94	-25.62	-30.21	-21.17	17.45	12.38	21.52
2011-08-04	20.12	18.69	21.68	-30.48	-37.29	-23.40	19.84	13.44	26.60

30 days fit total; GPP/ER/K600 all in reasonable ranges.

`predict_DO()` — downstream DO fit (first 10 rows)

solar.time	DO.obs.down	DO.mod.down
2011-07-31 00:00:00	6.95	6.90
2011-07-31 00:15:00	6.94	6.88
2011-07-31 00:30:00	6.95	6.86
2011-07-31 00:45:00	6.96	6.88
2011-07-31 01:00:00	6.96	6.85
2011-07-31 01:15:00	6.96	6.85
2011-07-31 01:30:00	6.97	6.86
2011-07-31 01:45:00	6.97	6.85
2011-07-31 02:00:00	6.98	6.87
2011-07-31 02:15:00	6.99	6.88

Modeled DO tracks observed within ~0.05-0.1 mg/L.

aappling-usgs

Amazing body of work here! I'm not quite done reviewing, but submitting a partial review tonight as promised.

When I come back, I'll finish up looking at the tail end of metab_2station.R and will look at the files mm_validate_data.R through NAMESPACE.

…bayes_2s Renames the two-station model constructor/class/file to metab_bayes_2s and the prep function to prepdata_bayes_2s (moved into metab_bayes_2s.R, marked internal), so get_params.metab_model's metab_<type>()-based constructor lookup works correctly. Clarifies the b2 token-matching comment in mm_parse_name.R. Spells out \"variable flow\" in the VFTS acronym expansion everywhere it appears in source docs (metab_bayes_2s.R, mm_name.R, data.R, two_station_example.R).

Moves the travel.time bounds check into mm_validate_data() (gated on travel.time's presence, so one-station data is unaffected). Adds mm_validate_data_2station() for the lead-in-coverage check, called early in metab_bayes_2s() before any data prep. Removes the now-unreachable duplicate lead-in check from prepdata_bayes_2s(). Per team discussion, the travel.time upper bound is tightened from <1 day to <=8 hours: beyond a units-mistake guard, this prevents the previous day's light conditions from influencing the following day's metabolism estimate. Error message and roxygen docs (metab_bayes_2s()'s Two-station data requirements section) updated to state both rationales. Resolves the scientific question raised by bishopia's review comment on lead-in duration/sunlight exposure.

Matches the convention used by all other default spec values, rather than computing fallback defaults inline in prepdata_bayes_2s().

Adds mm_time_by_date_matrix() and mm_check_dates_contiguous() as shared internal helpers, replacing near-duplicate inline logic in prepdata_bayes() and prepdata_bayes_2s(). Naming echoes Alison's review comment wording. Adds direct unit test coverage for prepdata_bayes() for the first time (previously only exercised indirectly via Stan-fitting integration tests). Fixes a real bug surfaced while extracting: prepdata_bayes()'s contiguous- sort check was missing both an isTRUE() wrapper and an is.list() guard present in prepdata_bayes_2s()'s version, causing it to fail with an opaque 'invalid argument type' error instead of its intended message whenever it should have fired. Both now share the corrected check.

Replaces near-identical detectCores()/fallback/min logic in runstan_bayes() and metab_bayes_2s() with a shared internal helper. Preserves runstan_bayes()'s original message. Adds docs/tests for mm_determine_cores()

metab_bayes_2s() previously treated a non-converged/failed Stan run (mode==2L) as a hard error, aborting post-processing via stop(). Now matches runstan_bayes()'s existing warn-and-continue pattern: emits a warning with the same diagnostic content, skips formatting/joining (which assume a successful fit), and lets the existing NA-fallback path (real dates, NA daily estimates, inst=NULL) handle the rest - unchanged logic, now reached via a caught warning instead of a caught error. The diagnostic now lands in fit\$warnings instead of fit\$errors, consistent with one-station's convention that a non-converged run is a recoverable, expected outcome rather than an internal failure. Checked for shareable overlap with runstan_bayes()'s equivalent handling before committing (per Alison's preference against duplication) - found differences that make extraction tricky at this point (before a larger refactor). Logged these for a future consolidation (when 2 station allows for separate fitting per day).

DO.obs.up/DO.sat.up/DO.obs.down/DO.sat.down docs: 'mg O2 / L' -> 'mgO2 L^-1', matching unitted-style notation used elsewhere.

@PARAM

@PARAM data_all previously claimed 'just one estimation-day' - actually receives the full filtered multi-day dataset when specs\$split_dates==FALSE. Fixed to accurately describe the bayes_1ply()/bayes_allply() split.

mdodrill-usgs · 2026-08-03T18:50:51Z

Ready for another look - summary of what's landed since your last review, grouped roughly by theme (individual replies have the full detail):

Naming/structure (a6eb23b): metab_2station → metab_bayes_2s, mm_ts_prep_data → prepdata_bayes_2s, moved into metab_bayes_2s.R, marked internal.

Validation (2074de0): consolidated into mm_validate_data_2station(); travel.time bound tightened to ≤8h per team discussion (light-exposure rationale, not just a units guard).

Shared logic extracted: matrix-pivot + contiguous-sort (5f3030e, fixed a bug along the way), core-count logic (4103b71). Stan-calling machinery itself not shared yet - current usage patterns genuinely differ (two-station always fits jointly, one station doesn't), but that's tied to the joint-fit constraint specifically, not a permanent stance. Once #475's per-day fitting work lands, sharing this becomes the natural next step, noted in that issue.

Behavior alignment (13b8c8f): failed Stan runs now warn-and-continue, matching one-station, instead of erroring.

Docs/cleanup: spec defaults centralized (8841358), units notation + stale bayes_allply() doc fixed (8a81e49).

Deferred to issue #475: day-filtering/gap tolerance, upstream covariates, per-day fitting, and a few smaller related questions - scoped out but decided this needs its own follow-up pass rather than partial implementation here.

All individual comment threads replied to. Most of them have also been marked resolved. Let me know if anything needs another pass.

mdodrill-usgs · 2026-08-20T20:02:19Z

Hi @aappling-usgs, just bumping this to the top. The PR is ready for another round of review - let me know if you think anything needs more attention.

aappling-usgs

^ GitHub congratulates me on reading to the end of your PR - I did it! More importantly, you did it! Huge and exciting body of changes here!

I did add a few more comments, but many fewer than before and some of them are thoughts for future PRs anyway. I'm approving now on the assumption that you'll think through these new comments and decide whether to resolve them here or in a future PR.

Mirrors the attach.units=deprecated() pattern already used in predict_metab.metab_bayes(): warns via unitted_deprecate_warn() when the argument is supplied, and attaches units via get_units(mm_data()) on request. mm_data() was missing unit entries for the bare K600/K600.lower/K600.upper column names (only K600.daily/.lower/.upper existed), which predict_metab's output actually uses -- added those (d^-1, matching the existing K600.daily convention) so the units lookup resolves. Addresses PR review comments on metab_bayes_2s.R:419.

Trimmed maintainer-only rationale out of rendered docs, relocated to code comments for future maintainers: - get_params.metab_bayes_2s(): moved the explanation of why this doesn't delegate via NextMethod() out of roxygen and into a comment above uncertainty <- match.arg(uncertainty). - mm_validate_data_2station(): cut a near-duplicate restatement of the lead-in-coverage requirement (already documented on metab_bayes_2s()'s own @section); now points to that as the source of truth and keeps only the mechanism note. Also includes streamMetabolizer.Rd's newly-rendered \seealso/\author block, surfaced by devtools::document() after fixing an unrelated pre-existing @doctype deprecation -- real package metadata that was silently missing from prior renders, not new content.

mdodrill-usgs · 2026-08-28T17:42:21Z

All review comments addressed — summary for anyone landing here later:

metab_bayes_2s.R: added attach.units deprecation support to predict_metab.metab_bayes_2s() (matching one-station's pattern); confirmed it can't fully delegate to predict_metab.metab_bayes() since that function doesn't handle K600 at all.
mm_validate_data.R: confirmed the travel.time hard bound stays as-is for now; the fraction-of-previous-day-light idea is captured as a future-dev item. Clarified the mm_validate_data()/mm_validate_data_2station() split and rationale.
specs.R: these threads surfaced a genuine unit mismatch in the two-station K600 prior (wrong physical scale relative to what the Stan model's math expects). Scoped this out (planned after feat/two-station-data-prep).
specs.R: agreed the bayes/bayes_2s switch() cases should converge once the generator work lands.

Also ran a doc-trim pass on this branch's own roxygen additions before closing out.

mdodrill-usgs added 7 commits July 9, 2026 13:17

Fix R CMD CHECK NOTEs in metab_2station.R

c76e7d6

- Add utils::globalVariables for ., metab_50pct, DO.mod.down - Add @importFrom utils modifyList - Replace require(rstan) with requireNamespace() pattern

aappling-usgs mentioned this pull request Jul 20, 2026

Provide helper functions for two-station methods #442

Open

aappling-usgs reviewed Jul 20, 2026

View reviewed changes

aappling-usgs mentioned this pull request Jul 21, 2026

Modify .stan generation to include two station option #441

Open

mdodrill-usgs mentioned this pull request Jul 21, 2026

Add alias scheme for arguments in R/mm_name.R #472

Open

mdodrill-usgs added 4 commits July 27, 2026 12:49

Move K600_lnorm_meanlog/sdlog defaults into specs()

8841358

Matches the convention used by all other default spec values, rather than computing fallback defaults inline in prepdata_bayes_2s().

mdodrill-usgs mentioned this pull request Jul 29, 2026

prepdata_bayes()'s "modal timestep" check doesn't do what its name/comment claim #474

Open

Extract shared core-count logic into mm_determine_cores()

4103b71

Replaces near-identical detectCores()/fallback/min logic in runstan_bayes() and metab_bayes_2s() with a shared internal helper. Preserves runstan_bayes()'s original message. Adds docs/tests for mm_determine_cores()

mdodrill-usgs mentioned this pull request Jul 31, 2026

Add test dataset for two-station #459

Open

mdodrill-usgs mentioned this pull request Aug 1, 2026

Two-station data prep: support real deployment data (gaps, multi-deployment, upstream covariates, per-day fitting) #475

Open

mdodrill-usgs added 2 commits August 2, 2026 09:43

Fix units notation to unitted-style

d270bfe

DO.obs.up/DO.sat.up/DO.obs.down/DO.sat.down docs: 'mg O2 / L' -> 'mgO2 L^-1', matching unitted-style notation used elsewhere.

Fix stale bayes_allply() data_all doc

8a81e49

@PARAM data_all previously claimed 'just one estimation-day' - actually receives the full filtered multi-day dataset when specs\$split_dates==FALSE. Fixed to accurately describe the bayes_1ply()/bayes_allply() split.

mdodrill-usgs requested a review from aappling-usgs August 3, 2026 18:52

mdodrill-usgs added the CDI26 label Aug 11, 2026

aappling-usgs approved these changes Aug 26, 2026

View reviewed changes

Comment thread R/metab_bayes_2s.R

Comment thread R/metab_bayes_2s.R

Comment thread R/mm_validate_data.R

Comment thread R/mm_validate_data.R

Comment thread R/specs.R

Comment thread R/specs.R

Comment thread R/specs.R

Comment thread tests/testthat/test-metab_bayes.R

mdodrill-usgs added 2 commits August 28, 2026 09:30

mdodrill-usgs merged commit eeb2a6c into dev/two-station Aug 28, 2026

mdodrill-usgs deleted the fix/vfts-stan-syntax branch August 28, 2026 22:47

Uh oh!

Conversation

mdodrill-usgs commented Jul 13, 2026 • edited by aappling-usgs Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Changes

Tests

Known gaps (deferred, noted in code)

References

Uh oh!

mdodrill-usgs commented Jul 13, 2026

Demo

Dataset

Model fit

predict_metab() — daily estimates (first 5 days)

predict_DO() — downstream DO fit (first 10 rows)

Uh oh!

aappling-usgs left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mdodrill-usgs commented Aug 3, 2026

Uh oh!

mdodrill-usgs commented Aug 20, 2026

Uh oh!

aappling-usgs left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mdodrill-usgs commented Aug 28, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

mdodrill-usgs commented Jul 13, 2026 •

edited by aappling-usgs

Loading

`predict_metab()` — daily estimates (first 5 days)

`predict_DO()` — downstream DO fit (first 10 rows)

aappling-usgs left a comment •

edited

Loading