From dce238b38bad11eb8d9693e2e315d07f0715dc58 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Thu, 9 Jul 2026 13:17:20 -0700 Subject: [PATCH 01/17] Add static two-station Stan model 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 --- inst/models/b2_np_oi_tr_plrckm.stan | 43 +++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 inst/models/b2_np_oi_tr_plrckm.stan diff --git a/inst/models/b2_np_oi_tr_plrckm.stan b/inst/models/b2_np_oi_tr_plrckm.stan new file mode 100644 index 0000000..a39442f --- /dev/null +++ b/inst/models/b2_np_oi_tr_plrckm.stan @@ -0,0 +1,43 @@ +// b2_np_oi_tr_plrckm.stan + +data { + int n_obs; // number of total do observations + int n_days; // number of days + array[n_obs] vector[n_days] DO_obs_up; + array[n_obs] vector[n_days] DO_sat_up; + array[n_obs] vector[n_days] DO_obs_down; + array[n_obs] vector[n_days] DO_sat_down; + array[n_obs] vector[n_days] light; + array[n_obs] vector[n_days] depth; + array[n_obs] vector[n_days] temp_water; + array[n_obs] vector[n_days] travel_time; + real K600_lnorm_meanlog; + real K600_lnorm_sdlog; +} + +parameters { + vector[n_days] GPP; + vector[n_days] ER; + vector[n_days] k600; + real sigma; +} + +transformed parameters { + array[n_obs] vector[n_days] metab; + array[n_obs] vector[n_days] KO2; + for (i in 1:n_obs){ + for (t in 1:n_days){ + KO2[i,t] = (k600[t] / depth[i,t]) / ((600 / (1800.6 - (temp_water[i,t] * 120.1) + (3.7818 * temp_water[i,t]^2) - (0.047608 * temp_water[i,t]^3)))^-0.5); + } + metab[i] = (DO_obs_up[i] + GPP .* light[i] ./ depth[i] + ER .* travel_time[i] ./ depth[i] + (KO2[i] .* travel_time[i] .* (DO_sat_up[i] - DO_obs_up[i] + DO_sat_down[i]) / 2)) ./ (1 + (KO2[i] .* travel_time[i]) / 2); + } +} + +model { + for (i in 1:n_obs){ + DO_obs_down[i] ~ normal(metab[i], sigma); + } + for (i in 1:n_days){ + k600[i] ~ lognormal(K600_lnorm_meanlog, K600_lnorm_sdlog); + } +} From b99fa34907adc080e240bc68980d01761db2bc06 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Thu, 9 Jul 2026 15:13:16 -0700 Subject: [PATCH 02/17] Add two-station data spec, validator stub, and constructor stub - 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 --- R/metab_2station.R | 67 ++++++++++++++++++++++++++++ R/mm_data.R | 23 ++++++++++ R/mm_validate_data.R | 5 ++- tests/testthat/test-metab_2station.R | 53 ++++++++++++++++++++++ 4 files changed, 147 insertions(+), 1 deletion(-) create mode 100644 R/metab_2station.R create mode 100644 tests/testthat/test-metab_2station.R diff --git a/R/metab_2station.R b/R/metab_2station.R new file mode 100644 index 0000000..cba7278 --- /dev/null +++ b/R/metab_2station.R @@ -0,0 +1,67 @@ +#' Two-station Bayesian metabolism model fitting function (stub) +#' +#' Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to +#' estimate GPP and ER from paired upstream and downstream DO, temperature, +#' light, and travel-time data. This function is currently a stub: it +#' validates the \code{data} argument and enforces two-station-specific data +#' requirements, but does not yet fit a model. See \code{\link{mm_name}} to +#' choose a Bayesian model and \code{\link{specs}} for relevant options for +#' the \code{specs} argument. +#' +#' @inheritParams metab +#' @return Not yet implemented; currently always errors after data validation. +#' +#' @section Two-station data requirements: In addition to the checks +#' performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the +#' reach travel time between the upstream and downstream stations, in days) +#' must be strictly positive and less than 1 day (values >= 1 usually +#' indicate travel time was supplied in the wrong units). There must also be +#' enough lead-in observations of upstream DO before the first row of +#' \code{data} to cover the longest travel time in the dataset, given the +#' (median) timestep of \code{data$solar.time}. +#' +#' @export +#' @family metab_model +metab_2station <- function( + specs=specs(mm_name('bayes_2station')), + data=mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, + light, depth, temp.water, travel.time), + data_daily=mm_data(date, optional='all'), + info=NULL +) { + + # Check data for correct column names & units + dat_list <- mm_validate_data(data, data_daily, 'metab_2station') + + travel_time <- v(dat_list$data$travel.time) + solar_time <- v(dat_list$data$solar.time) + + # a. travel.time must be strictly positive + if(any(travel_time <= 0)) { + stop('travel.time must be > 0') + } + + # b. travel.time is expected in days, so any value >= 1 almost certainly + # reflects a units mistake (e.g., minutes or hours rather than days) + if(any(travel_time >= 1)) { + stop('travel.time must be < 1 day; values >= 1 suggest incorrect units (expected days)') + } + + # c. there must be enough lead-in rows of upstream DO before the first + # modeled row to cover the longest travel time in the dataset. timestep_days + # is the median observation interval, in days; max_lag is the number of + # timesteps by which upstream data must lead downstream predictions. The + # first max_lag rows of data serve only as lead-in and cannot themselves be + # modeled, so at least max_lag + 1 rows are required overall. + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + max_lag <- max(round(travel_time / timestep_days)) + if(nrow(dat_list$data) <= max_lag) { + lead_in_needed <- max_lag - nrow(dat_list$data) + 1 + stop(paste0( + 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', + max_lag, ' timestep(s), but only ', nrow(dat_list$data), ' row(s) were supplied; ', + 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) + } + + stop('metab_2station not yet implemented') +} diff --git a/R/mm_data.R b/R/mm_data.R index d9ce239..e6009c4 100644 --- a/R/mm_data.R +++ b/R/mm_data.R @@ -24,6 +24,24 @@ #' equilibrium saturation \eqn{mg O[2] L^{-1}}{mg O2 / L}. Calculate using #' \link{calc_DO_sat}} #' +#' \item{ \code{DO.obs.up} dissolved oxygen concentration observations at the +#' upstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{DO.sat.up} dissolved oxygen concentrations at equilibrium +#' saturation at the upstream station of a two-station reach, \eqn{mg O[2] +#' L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{DO.obs.down} dissolved oxygen concentration observations at +#' the downstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 +#' / L}} +#' +#' \item{ \code{DO.sat.down} dissolved oxygen concentrations at equilibrium +#' saturation at the downstream station of a two-station reach, \eqn{mg O[2] +#' L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{travel.time} reach travel time between the upstream and +#' downstream stations of a two-station reach, in days, \eqn{d}{d}} +#' #' \item{ \code{depth} stream depth, \eqn{m}{m}}. #' #' \item{ \code{temp.water} water temperature, \eqn{degC}}. @@ -97,9 +115,14 @@ mm_data <- function(..., optional='none') { solar.time = u(as.POSIXct("2050-03-14 15:10:00", tz="UTC"), NA), DO.obs = u(10.1,"mgO2 L^-1"), DO.sat = u(14.2,"mgO2 L^-1"), + DO.obs.up = u(10.1,"mgO2 L^-1"), + DO.sat.up = u(14.2,"mgO2 L^-1"), + DO.obs.down = u(9.8,"mgO2 L^-1"), + DO.sat.down = u(14.0,"mgO2 L^-1"), depth = u(0.5,"m"), temp.water = u(21.8,"degC"), light = u(300.9,"umol m^-2 s^-1"), + travel.time = u(0.05,"d"), discharge = u(9,"m^3 s^-1"), velocity = u(2,"m s^-1"), date = u(as.Date("2050-03-14", tz="UTC"), NA), diff --git a/R/mm_validate_data.R b/R/mm_validate_data.R index c9fbea4..401b9bc 100644 --- a/R/mm_validate_data.R +++ b/R/mm_validate_data.R @@ -59,7 +59,10 @@ mm_validate_data <- function( # missing_cols was not among the data_tests or the metab_model data were # specified without a timestamp column if('na_times' %in% data_tests) { - timecol <- grep('date|time', names(dat), value=TRUE) + # match against the known timestamp column names rather than a + # substring grep for 'date'/'time', which would also match non- + # timestamp columns such as 'travel.time' + timecol <- intersect(c('solar.time','date'), names(dat)) if(length(timecol) != 1) stop("in ", data_type, " found ", length(timecol), " possible timestamp columns", call.=FALSE) na.times <- which(is.na(dat[[timecol]])) if(length(na.times) > 0) { diff --git a/tests/testthat/test-metab_2station.R b/tests/testthat/test-metab_2station.R new file mode 100644 index 0000000..ec546ad --- /dev/null +++ b/tests/testthat/test-metab_2station.R @@ -0,0 +1,53 @@ + +# Build a minimal, valid two-station data.frame. Defaults give a 5-minute +# timestep (0.0034722 days) and a 0.01-day travel time, so +# max_lag = round(0.01 / 0.0034722) = 3 timesteps of required upstream lead-in. +make_2station_data <- function(n=10, timestep_min=5, travel_time=0.01) { + data.frame( + solar.time = as.POSIXct("2050-06-01 00:00:00", tz="UTC") + + as.difftime((seq_len(n) - 1) * timestep_min, units="mins"), + DO.obs.up = rep(9, n), + DO.sat.up = rep(10, n), + DO.obs.down = rep(8.8, n), + DO.sat.down = rep(9.9, n), + light = rep(300, n), + depth = rep(0.5, n), + temp.water = rep(20, n), + travel.time = rep(travel_time, n) + ) +} + +test_that("mm_validate_data catches missing required columns", { + dat <- dplyr::select(make_2station_data(), -DO.obs.up) + expect_error(metab_2station(data=dat), "missing these columns") +}) + +test_that("travel.time <= 0 triggers an error", { + dat <- make_2station_data(travel_time=0) + expect_error(metab_2station(data=dat), "travel.time must be > 0") + + dat <- make_2station_data(travel_time=-0.01) + expect_error(metab_2station(data=dat), "travel.time must be > 0") +}) + +test_that("travel.time >= 1 triggers an error with a units hint", { + dat <- make_2station_data(travel_time=1) + expect_error(metab_2station(data=dat), "travel.time must be < 1 day.*incorrect units") + + dat <- make_2station_data(travel_time=1.5) + expect_error(metab_2station(data=dat), "travel.time must be < 1 day.*incorrect units") +}) + +test_that("insufficient lead-in data triggers an error", { + # 2 rows but max_lag=3 timesteps of upstream lead-in are needed + dat <- make_2station_data(n=2) + expect_error(metab_2station(data=dat), "insufficient lead-in data") +}) + +test_that("a minimal valid data.frame passes all data-format checks", { + dat <- make_2station_data(n=10) + # metab_2station is a stub, so a fully valid data.frame should make it all + # the way through the column/units/travel.time/lead-in checks and fail only + # on the final "not yet implemented" stop -- not on any validation check. + expect_error(metab_2station(data=dat), "metab_2station not yet implemented") +}) From d8a1d04b9f1f5bf21fbfdf83fb9d17743893dd16 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Thu, 9 Jul 2026 15:52:39 -0700 Subject: [PATCH 03/17] Add mm_ts_prep_data() with upstream lag shift - 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 --- R/mm_ts_prep_data.R | 125 +++++++++++++++++++++++++++ tests/testthat/test-metab_2station.R | 99 +++++++++++++++++++++ 2 files changed, 224 insertions(+) create mode 100644 R/mm_ts_prep_data.R diff --git a/R/mm_ts_prep_data.R b/R/mm_ts_prep_data.R new file mode 100644 index 0000000..5414261 --- /dev/null +++ b/R/mm_ts_prep_data.R @@ -0,0 +1,125 @@ +#' Reshape long-format two-station data into the list expected by the +#' two-station Stan model +#' +#' Time-shifts the upstream DO series to match the travel time between +#' stations, then pivots the result into the \code{n_obs x n_days} matrices +#' expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} +#' (see \code{\link{metab_2station}}). +#' +#' The upstream observation that "matches" a given downstream observation at +#' row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / +#' timestep_days)} timesteps earlier, where \code{timestep_days} is the +#' median timestep of \code{data$solar.time}. This must be computed the same +#' way as in \code{\link{metab_2station}}'s lead-in check, so that the +#' \code{max(lag)} computed here always agrees with the lead-in requirement +#' already validated there. The first \code{max(lag)} rows of \code{data} are +#' lead-in rows: they supply upstream DO for the shift but are never +#' themselves treated as modeled (downstream) observations. +#' +#' @param data data.frame as validated by \code{\link{mm_validate_data}} for +#' \code{\link{metab_2station}}: must contain \code{solar.time}, +#' \code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, +#' \code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, +#' \code{travel.time}, sorted ascending by \code{solar.time}, and must +#' include the lead-in rows required to cover the longest travel time (see +#' \code{\link{metab_2station}}). +#' @param specs optional list of model specs. If it contains +#' \code{K600_lnorm_meanlog} and/or \code{K600_lnorm_sdlog}, those values +#' are used; otherwise placeholder defaults of \code{log(3.48)} and +#' \code{0.5} are used. +#' @return a named list with all variables in the Stan model's data block: +#' \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, +#' \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, +#' \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} +#' matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} +#' @importFrom unitted v +#' @export +mm_ts_prep_data <- function(data, specs=NULL) { + + # strip units; Stan cannot handle unitted vectors/matrices + data <- v(data) + + # timestep_days must match metab_2station()'s lead-in check exactly (median + # timestep, in days), so that max_lag here agrees with what was validated + # there + timestep_days <- stats::median(as.numeric(diff(data$solar.time), units='days')) + + # lag, in timesteps, that the upstream series must be shifted by to line up + # with each row's downstream observation + lag <- round(data$travel.time / timestep_days) + max_lag <- max(lag) + n_total <- nrow(data) + if(n_total <= max_lag) { + stop( + 'not enough lead-in rows to cover the longest travel.time (', max_lag, ' timesteps ', + 'implied, but only ', n_total, ' rows supplied); this should already have been caught ', + 'by metab_2station()') + } + + # the first max_lag rows are lead-in only (upstream DO used for the shift, + # but never modeled themselves). every row i > max_lag is guaranteed to + # have a valid shift target (i - lag[i] >= 1) because lag[i] <= max_lag + keep <- seq.int(max_lag + 1, n_total) + shift_idx <- keep - lag[keep] + if(any(shift_idx < 1)) { + # should be unreachable given max_lag's definition; guards against + # programming errors rather than expected user input + stop('internal error: upstream shift index falls before the first row of data') + } + + modeled <- data.frame( + solar.time = data$solar.time[keep], + DO_obs_up = data$DO.obs.up[shift_idx], + DO_sat_up = data$DO.sat.up[shift_idx], + DO_obs_down = data$DO.obs.down[keep], + DO_sat_down = data$DO.sat.down[keep], + light = data$light[keep], + depth = data$depth[keep], + temp_water = data$temp.water[keep], + travel_time = data$travel.time[keep] + ) + + # pivot into n_obs x n_days matrices, one column per unique date, following + # the same time_by_date_matrix approach used for the 1-station Stan models + # (see prepdata_bayes() in metab_bayes.R) + date_vec <- as.character(as.Date(modeled$solar.time)) + date_table <- table(date_vec) + n_days <- length(date_table) + n_obs_per_day <- unique(unname(date_table)) + if(length(n_obs_per_day) > 1) { + stop( + 'dates have differing numbers of modeled rows after lead-in removal; ', + 'observations cannot be combined into a matrix: ', + paste(sprintf('%s (%d rows)', names(date_table), date_table), collapse=', ')) + } + n_obs <- n_obs_per_day + + to_matrix <- function(vec) matrix(vec, nrow=n_obs, ncol=n_days, byrow=FALSE) + + # confirm each date occupies a contiguous block of rows, i.e., that data + # was sorted by solar.time; otherwise the matrix pivot below would silently + # scramble which rows belong to which date + date_mat <- to_matrix(date_vec) + unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) + if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { + stop('data must be sorted by solar.time so that each date occupies a contiguous block of rows') + } + + K600_lnorm_meanlog <- if(!is.null(specs$K600_lnorm_meanlog)) specs$K600_lnorm_meanlog else log(3.48) + K600_lnorm_sdlog <- if(!is.null(specs$K600_lnorm_sdlog)) specs$K600_lnorm_sdlog else 0.5 + + list( + n_obs = n_obs, + n_days = n_days, + DO_obs_up = to_matrix(modeled$DO_obs_up), + DO_sat_up = to_matrix(modeled$DO_sat_up), + DO_obs_down = to_matrix(modeled$DO_obs_down), + DO_sat_down = to_matrix(modeled$DO_sat_down), + light = to_matrix(modeled$light), + depth = to_matrix(modeled$depth), + temp_water = to_matrix(modeled$temp_water), + travel_time = to_matrix(modeled$travel_time), + K600_lnorm_meanlog = K600_lnorm_meanlog, + K600_lnorm_sdlog = K600_lnorm_sdlog + ) +} diff --git a/tests/testthat/test-metab_2station.R b/tests/testthat/test-metab_2station.R index ec546ad..1d1c281 100644 --- a/tests/testthat/test-metab_2station.R +++ b/tests/testthat/test-metab_2station.R @@ -51,3 +51,102 @@ test_that("a minimal valid data.frame passes all data-format checks", { # on the final "not yet implemented" stop -- not on any validation check. expect_error(metab_2station(data=dat), "metab_2station not yet implemented") }) + + +# mm_ts_prep_data() ----------------------------------------------------- + +# Build a two-day, unit-labeled data.frame with a known, traceable +# DO.obs.up/DO.sat.up series (sequential integers) so the shift can be +# checked by exact value, plus a leading lead-in block. 5-minute timestep +# (0.0034722 days) and 0.01-day travel.time give +# max_lag = round(0.01 / 0.0034722) = 3 lead-in timesteps. +# Day 1 = 10 rows (3 lead-in + 7 modeled), Day 2 = 7 rows (all modeled), so +# both modeled days end up with n_obs = 7 rows. +make_ts_data <- function(n_leadin=3, n_day1=10, n_day2=7, travel_time=0.01, unitted=FALSE) { + n_total <- n_day1 + n_day2 + solar.time <- c( + as.POSIXct("2050-06-01 00:00:00", tz="UTC") + as.difftime((seq_len(n_day1) - 1) * 5, units="mins"), + as.POSIXct("2050-06-02 00:00:00", tz="UTC") + as.difftime((seq_len(n_day2) - 1) * 5, units="mins")) + dat <- data.frame( + solar.time = solar.time, + DO.obs.up = seq_len(n_total), # traceable: value == original row index + DO.sat.up = seq_len(n_total) + 100, # traceable, offset so it's distinguishable from DO.obs.up + DO.obs.down = seq_len(n_total) + 1000, # traceable, offset so it's distinguishable from up/sat values + DO.sat.down = rep(9.9, n_total), + light = rep(300, n_total), + depth = rep(0.5, n_total), + temp.water = rep(20, n_total), + travel.time = rep(travel_time, n_total) + ) + if(unitted) { + units_template <- get_units(mm_data( + solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, depth, temp.water, travel.time)) + dat <- u(dat, unname(units_template[names(dat)])) + } + dat +} + +test_that("upstream DO is shifted by the correct lag", { + dat <- make_ts_data() + out <- mm_ts_prep_data(dat) + + # max_lag=3, so modeled row i (original index i) uses upstream data from + # original row (i - 3). day 1's 7 modeled rows are original rows 4:10, so + # they pick up DO.obs.up from original rows 1:7; day 2's 7 modeled rows are + # original rows 11:17, picking up DO.obs.up from original rows 8:14. + expect_equal(out$DO_obs_up[,1], as.numeric(1:7)) + expect_equal(out$DO_obs_up[,2], as.numeric(8:14)) + expect_equal(out$DO_sat_up[,1], as.numeric(1:7) + 100) + expect_equal(out$DO_sat_up[,2], as.numeric(8:14) + 100) +}) + +test_that("lead-in rows are excluded from the output matrices", { + dat <- make_ts_data(n_leadin=3, n_day1=10, n_day2=7) + out <- mm_ts_prep_data(dat) + + # 17 total rows in, 3 are lead-in-only, so 14 modeled rows should remain + expect_equal(out$n_obs * out$n_days, nrow(dat) - 3) + # none of the lead-in DO.obs.up values (1, 2, 3) should appear as a + # DOWNSTREAM-paired value, i.e., the first modeled column should start at + # the shifted value 1, not 1:3 appearing as downstream/lead-in rows + expect_false(any(dat$DO.obs.down[1:3] %in% unlist(out$DO_obs_down))) +}) + +test_that("output matrices have n_obs x n_days dimensions", { + dat <- make_ts_data() + out <- mm_ts_prep_data(dat) + + expect_equal(out$n_obs, 7) + expect_equal(out$n_days, 2) + for(varname in c('DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down','light','depth','temp_water','travel_time')) { + expect_equal(dim(out[[varname]]), c(7, 2), info=varname) + } +}) + +test_that("all required Stan data block variables are present", { + dat <- make_ts_data() + out <- mm_ts_prep_data(dat) + + expected_names <- c( + 'n_obs','n_days','DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down', + 'light','depth','temp_water','travel_time','K600_lnorm_meanlog','K600_lnorm_sdlog') + expect_true(all(expected_names %in% names(out))) + + # placeholder K600 lognormal priors, per PR D-3 + expect_equal(out$K600_lnorm_meanlog, log(3.48)) + expect_equal(out$K600_lnorm_sdlog, 0.5) +}) + +test_that("units are stripped from all numeric outputs", { + dat <- make_ts_data(unitted=TRUE) + expect_true(is.unitted(dat)) + + out <- mm_ts_prep_data(dat) + for(varname in c('DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down','light','depth','temp_water','travel_time')) { + expect_false(is.unitted(out[[varname]]), info=varname) + } + expect_false(is.unitted(out$K600_lnorm_meanlog)) + expect_false(is.unitted(out$K600_lnorm_sdlog)) + expect_false(is.unitted(out$n_obs)) + expect_false(is.unitted(out$n_days)) +}) From 8246c917ac685febfc0fc49ed0a00f75447e5061 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Fri, 10 Jul 2026 12:53:33 -0700 Subject: [PATCH 04/17] Add predict_DO.metab_2station stub and fix mm_parse_name for b2_ prefix - 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 --- DESCRIPTION | 4 +- NAMESPACE | 7 ++- R/metab_2station.R | 19 ++++++++ R/mm_parse_name.R | 6 ++- man/metab_2station.Rd | 68 ++++++++++++++++++++++++++++ man/mm_ts_prep_data.Rd | 47 +++++++++++++++++++ tests/testthat/test-metab_2station.R | 40 ++++++++++++++++ 7 files changed, 188 insertions(+), 3 deletions(-) create mode 100644 man/metab_2station.Rd create mode 100644 man/mm_ts_prep_data.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 4af05e3..acf9b7d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -88,6 +88,7 @@ Collate: 'lookup_timezone.R' 'lookup_usgs_elevation.R' 'metab.R' + 'metab_2station.R' 'metab_model_interface.R' 'specs-class.R' 'metab_model-class.R' @@ -120,6 +121,7 @@ Collate: 'mm_predict_DO_1ply.R' 'mm_predict_metab_1ply.R' 'mm_sd_to_ci.R' + 'mm_ts_prep_data.R' 'mm_validate_data.R' 'plot_DO_preds.R' 'plot_distribs.R' @@ -129,5 +131,5 @@ Collate: 'streamMetabolizer-deprecated.R' 'streamMetabolizer.R' 'zz_build_docs.R' -RoxygenNote: 7.1.1 +RoxygenNote: 7.3.3 Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index b161194..a19b9e5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ S3method(get_params,metab_model) S3method(get_params,metab_sim) S3method(get_specs,metab_model) S3method(get_version,metab_model) +S3method(predict_DO,metab_2station) S3method(predict_DO,metab_Kmodel) S3method(predict_DO,metab_model) S3method(predict_DO,metab_night) @@ -67,6 +68,7 @@ export(lookup_google_timezone) export(lookup_timezone) export(lookup_usgs_elevation) export(metab) +export(metab_2station) export(metab_Kmodel) export(metab_bayes) export(metab_inputs) @@ -82,6 +84,7 @@ export(mm_model_by_ply) export(mm_model_by_ply_prototype) export(mm_name) export(mm_parse_name) +export(mm_ts_prep_data) export(mm_valid_names) export(mm_validate_data) export(mm_validate_name) @@ -114,7 +117,6 @@ importFrom(LakeMetabolizer,sw.to.par.base) importFrom(graphics,abline) importFrom(graphics,plot) importFrom(graphics,points) -importFrom(lazyeval,lazy_dots) importFrom(lifecycle,deprecate_warn) importFrom(lifecycle,deprecated) importFrom(lifecycle,is_present) @@ -124,6 +126,9 @@ importFrom(lubridate,is.Date) importFrom(lubridate,is.POSIXct) importFrom(lubridate,tz) importFrom(lubridate,with_tz) +importFrom(rlang,as_name) +importFrom(rlang,enquos) +importFrom(rlang,quo_is_null) importFrom(stats,approx) importFrom(stats,approxfun) importFrom(stats,coef) diff --git a/R/metab_2station.R b/R/metab_2station.R index cba7278..36efca7 100644 --- a/R/metab_2station.R +++ b/R/metab_2station.R @@ -65,3 +65,22 @@ metab_2station <- function( stop('metab_2station not yet implemented') } + +#' @describeIn predict_DO Stub for two-station (VFTS) models. Not yet +#' functional: \code{metab_2station} does not yet produce a fitted Stan +#' model to predict from (see \code{\link{metab_2station}}), so this always +#' errors. Once implemented, this method will return a data.frame +#' with columns \code{solar.time}, \code{DO.obs.down} (the observed +#' downstream DO from the input data), and \code{DO.mod.down} (the +#' two-station Stan model's predicted downstream DO) -- unlike the +#' one-station \code{predict_DO} methods, which return \code{DO.obs}/ +#' \code{DO.mod}. +#' @export +predict_DO.metab_2station <- function(metab_model, date_start=NA, date_end=NA, ..., use_saved=TRUE) { + + # NOTE: R/plot_DO_preds.R and tests/testthat/helper-rmse_DO.R both + # hard-code the one-station DO.obs/DO.mod column names; they'll need to + # branch on (or be parameterized for) DO.obs.down/DO.mod.down once this + # method actually returns two-station predictions. + stop("predict_DO for two-station models not yet implemented — requires completion") +} diff --git a/R/mm_parse_name.R b/R/mm_parse_name.R index 43cd18a..1d8d824 100644 --- a/R/mm_parse_name.R +++ b/R/mm_parse_name.R @@ -37,7 +37,11 @@ mm_parse_name <- function(model_name, expand=FALSE) { # parse the name parsed <- strsplit(basename(model_name), "_|\\.") sapply(1:length(parsed), function(pnum) if(length(parsed[[pnum]]) <= 5) stop('missing one or more pieces in name: ', model_name[pnum])) - type <- unname(c(b='bayes', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) + # 'b2' (not just 'b') is the whole first token for two-station model files + # (e.g., 'b2_np_oi_tr_plrckm.stan'), since strsplit on "_|\\." above never + # splits within a token; an exact-match entry here is therefore sufficient + # and requires no change to the token-extraction logic itself + type <- unname(c(b='bayes', b2='bayes_2station', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) pool_K600 <- unname(c( np='none', Kn='normal', Kn0='normal_sdzero', Knx='normal_sdfixed', diff --git a/man/metab_2station.Rd b/man/metab_2station.Rd new file mode 100644 index 0000000..e7ba103 --- /dev/null +++ b/man/metab_2station.Rd @@ -0,0 +1,68 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metab_2station.R +\name{metab_2station} +\alias{metab_2station} +\title{Two-station Bayesian metabolism model fitting function (stub)} +\usage{ +metab_2station( + specs = specs(mm_name("bayes_2station")), + data = mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, + depth, temp.water, travel.time), + data_daily = mm_data(date, optional = "all"), + info = NULL +) +} +\arguments{ +\item{specs}{a list of model specifications and parameters for a model. +Although this may be specified manually (it's just a list), it is easier +and safer to use \code{\link{specs}} to generate the list, because the set +of required parameters and their defaults depends on the model given in the +\code{model_name} argument to \code{specs}. The help file for +\code{\link{specs}} lists the necessary parameters, describes them in +detail, and gives default values.} + +\item{data}{data.frame (not a tbl_df) of input data at the temporal +resolution of raw observations (unit-value). Columns must have the same +names, units, and format as the default. The solar.time column must also +have a timezone code ('tzone' attribute) of 'UTC'. See the +\strong{'Formatting \code{data}'} section below for a full description.} + +\item{data_daily}{data.frame containing inputs with a daily timestep. See the +\strong{'Formatting \code{data_daily}'} section below for a full +description.} + +\item{info}{any information, in any format, that you would like to store +within the metab_model object} +} +\value{ +Not yet implemented; currently always errors after data validation. +} +\description{ +Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to +estimate GPP and ER from paired upstream and downstream DO, temperature, +light, and travel-time data. This function is currently a stub: it +validates the \code{data} argument and enforces two-station-specific data +requirements, but does not yet fit a model. See \code{\link{mm_name}} to +choose a Bayesian model and \code{\link{specs}} for relevant options for +the \code{specs} argument. +} +\section{Two-station data requirements}{ + In addition to the checks + performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the + reach travel time between the upstream and downstream stations, in days) + must be strictly positive and less than 1 day (values >= 1 usually + indicate travel time was supplied in the wrong units). There must also be + enough lead-in observations of upstream DO before the first row of + \code{data} to cover the longest travel time in the dataset, given the + (median) timestep of \code{data$solar.time}. +} + +\seealso{ +Other metab_model: +\code{\link{metab_Kmodel}}, +\code{\link{metab_bayes}}, +\code{\link{metab_mle}}, +\code{\link{metab_night}}, +\code{\link{metab_sim}} +} +\concept{metab_model} diff --git a/man/mm_ts_prep_data.Rd b/man/mm_ts_prep_data.Rd new file mode 100644 index 0000000..c85f226 --- /dev/null +++ b/man/mm_ts_prep_data.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_ts_prep_data.R +\name{mm_ts_prep_data} +\alias{mm_ts_prep_data} +\title{Reshape long-format two-station data into the list expected by the +two-station Stan model} +\usage{ +mm_ts_prep_data(data, specs = NULL) +} +\arguments{ +\item{data}{data.frame as validated by \code{\link{mm_validate_data}} for +\code{\link{metab_2station}}: must contain \code{solar.time}, +\code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, +\code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, +\code{travel.time}, sorted ascending by \code{solar.time}, and must +include the lead-in rows required to cover the longest travel time (see +\code{\link{metab_2station}}).} + +\item{specs}{optional list of model specs. If it contains +\code{K600_lnorm_meanlog} and/or \code{K600_lnorm_sdlog}, those values +are used; otherwise placeholder defaults of \code{log(3.48)} and +\code{0.5} are used.} +} +\value{ +a named list with all variables in the Stan model's data block: + \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, + \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, + \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} + matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} +} +\description{ +Time-shifts the upstream DO series to match the travel time between +stations, then pivots the result into the \code{n_obs x n_days} matrices +expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} +(see \code{\link{metab_2station}}). +} +\details{ +The upstream observation that "matches" a given downstream observation at +row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / +timestep_days)} timesteps earlier, where \code{timestep_days} is the +median timestep of \code{data$solar.time}. This must be computed the same +way as in \code{\link{metab_2station}}'s lead-in check, so that the +\code{max(lag)} computed here always agrees with the lead-in requirement +already validated there. The first \code{max(lag)} rows of \code{data} are +lead-in rows: they supply upstream DO for the shift but are never +themselves treated as modeled (downstream) observations. +} diff --git a/tests/testthat/test-metab_2station.R b/tests/testthat/test-metab_2station.R index 1d1c281..aa376dd 100644 --- a/tests/testthat/test-metab_2station.R +++ b/tests/testthat/test-metab_2station.R @@ -150,3 +150,43 @@ test_that("units are stripped from all numeric outputs", { expect_false(is.unitted(out$n_obs)) expect_false(is.unitted(out$n_days)) }) + + +# mm_parse_name() for two-station models --------------------------------- + +test_that("mm_parse_name recognizes the b2_ prefix for two-station models", { + parsed <- mm_parse_name('b2_np_oi_tr_plrckm.stan') + + expect_equal(parsed$type, 'bayes_2station') + # the rest of the name is shared syntax with one-station bayes models and + # should parse the same way regardless of the b vs. b2 prefix + expect_equal(parsed$pool_K600, 'none') + expect_true(parsed$err_obs_iid) + expect_false(parsed$err_proc_acor) + expect_false(parsed$err_proc_iid) + expect_false(parsed$err_proc_GPP) + expect_equal(parsed$ode_method, 'trapezoid') + expect_equal(parsed$GPP_fun, 'linlight') + expect_equal(parsed$ER_fun, 'constant') + expect_equal(parsed$deficit_src, 'DO_mod') + expect_equal(parsed$engine, 'stan') + + # a one-station name with the same suffix should still parse as plain 'bayes' + expect_equal(mm_parse_name('b_np_oi_tr_plrckm.stan')$type, 'bayes') +}) + + +# predict_DO.metab_2station stub ----------------------------------------- + +test_that("predict_DO dispatches to the metab_2station stub and errors as expected", { + # metab_2station() is itself a stub (see tests above) and never returns a + # fitted model object, so there's no way yet to construct a real + # metab_2station instance. Mock the class attribute alone to exercise S3 + # dispatch of predict_DO() to predict_DO.metab_2station(). + mm <- structure(list(), class = c('metab_2station', 'metab_model')) + + expect_error( + predict_DO(mm), + "predict_DO for two-station models not yet implemented .* requires D-4 completion" + ) +}) From d8cb70e041b621df0171bc1c19bf2c98200f1eac Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Fri, 10 Jul 2026 13:17:20 -0700 Subject: [PATCH 05/17] Update two-station Stan model: parameter naming and priors - 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 --- inst/models/b2_np_oi_tr_plrckm.stan | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/inst/models/b2_np_oi_tr_plrckm.stan b/inst/models/b2_np_oi_tr_plrckm.stan index a39442f..bf2ce4b 100644 --- a/inst/models/b2_np_oi_tr_plrckm.stan +++ b/inst/models/b2_np_oi_tr_plrckm.stan @@ -11,14 +11,18 @@ data { array[n_obs] vector[n_days] depth; array[n_obs] vector[n_days] temp_water; array[n_obs] vector[n_days] travel_time; + real GPP_daily_mu; + real GPP_daily_sigma; + real ER_daily_mu; + real ER_daily_sigma; real K600_lnorm_meanlog; real K600_lnorm_sdlog; } parameters { - vector[n_days] GPP; - vector[n_days] ER; - vector[n_days] k600; + vector[n_days] GPP_daily; + vector[n_days] ER_daily; + vector[n_days] K600_daily; real sigma; } @@ -27,17 +31,19 @@ transformed parameters { array[n_obs] vector[n_days] KO2; for (i in 1:n_obs){ for (t in 1:n_days){ - KO2[i,t] = (k600[t] / depth[i,t]) / ((600 / (1800.6 - (temp_water[i,t] * 120.1) + (3.7818 * temp_water[i,t]^2) - (0.047608 * temp_water[i,t]^3)))^-0.5); + KO2[i,t] = (K600_daily[t] / depth[i,t]) / ((600 / (1800.6 - (temp_water[i,t] * 120.1) + (3.7818 * temp_water[i,t]^2) - (0.047608 * temp_water[i,t]^3)))^-0.5); } - metab[i] = (DO_obs_up[i] + GPP .* light[i] ./ depth[i] + ER .* travel_time[i] ./ depth[i] + (KO2[i] .* travel_time[i] .* (DO_sat_up[i] - DO_obs_up[i] + DO_sat_down[i]) / 2)) ./ (1 + (KO2[i] .* travel_time[i]) / 2); + metab[i] = (DO_obs_up[i] + GPP_daily .* light[i] ./ depth[i] + ER_daily .* travel_time[i] ./ depth[i] + (KO2[i] .* travel_time[i] .* (DO_sat_up[i] - DO_obs_up[i] + DO_sat_down[i]) / 2)) ./ (1 + (KO2[i] .* travel_time[i]) / 2); } } model { + GPP_daily ~ normal(GPP_daily_mu, GPP_daily_sigma); + ER_daily ~ normal(ER_daily_mu, ER_daily_sigma); + for (i in 1:n_days){ + K600_daily[i] ~ lognormal(K600_lnorm_meanlog, K600_lnorm_sdlog); + } for (i in 1:n_obs){ DO_obs_down[i] ~ normal(metab[i], sigma); } - for (i in 1:n_days){ - k600[i] ~ lognormal(K600_lnorm_meanlog, K600_lnorm_sdlog); - } } From 033e9d673edfb52f7ff070f36d2f54e506652cf3 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 13 Jul 2026 11:21:12 -0700 Subject: [PATCH 06/17] Implement metab_2station() with specs, dispatch, example data, and tests - 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 --- DESCRIPTION | 11 +- NAMESPACE | 3 + R/data.R | 44 +++ R/metab.R | 11 +- R/metab_2station.R | 382 +++++++++++++++++++++++---- R/mm_name.R | 33 ++- R/mm_valid_names.R | 7 +- R/specs.R | 71 ++++- data-raw/two_station_example.R | 89 +++++++ data/two_station_example.rda | Bin 0 -> 17629 bytes man/metab_2station-class.Rd | 26 ++ man/two_station_example.Rd | 55 ++++ tests/testthat/test-metab_2station.R | 87 ++++-- 13 files changed, 739 insertions(+), 80 deletions(-) create mode 100644 R/data.R create mode 100644 data-raw/two_station_example.R create mode 100644 data/two_station_example.rda create mode 100644 man/metab_2station-class.Rd create mode 100644 man/two_station_example.Rd diff --git a/DESCRIPTION b/DESCRIPTION index acf9b7d..d3852d2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -39,9 +39,12 @@ Imports: lubridate, magrittr, methods, + parallel, + stats, tibble (>= 1.1.0), tidyr, - unitted (>= 0.2.8) + unitted (>= 0.2.8), + utils Suggests: chron, devtools, @@ -80,6 +83,7 @@ Collate: 'create_calc_DO.R' 'create_calc_NLL.R' 'create_calc_dDOdt.R' + 'data.R' 'data_metab.R' 'deprecated.R' 'load_french_creek.R' @@ -88,12 +92,12 @@ Collate: 'lookup_timezone.R' 'lookup_usgs_elevation.R' 'metab.R' - 'metab_2station.R' 'metab_model_interface.R' 'specs-class.R' 'metab_model-class.R' - 'metab_Kmodel.R' 'metab_bayes.R' + 'metab_2station.R' + 'metab_Kmodel.R' 'metab_inputs.R' 'metab_mle.R' 'metab_model.get_param_names.R' @@ -133,3 +137,4 @@ Collate: 'zz_build_docs.R' RoxygenNote: 7.3.3 Encoding: UTF-8 +LazyData: true diff --git a/NAMESPACE b/NAMESPACE index a19b9e5..fc43e76 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,7 @@ S3method(get_mcmc,metab_bayes) S3method(get_mcmc_data,metab_bayes) S3method(get_param_names,character) S3method(get_param_names,metab_model) +S3method(get_params,metab_2station) S3method(get_params,metab_Kmodel) S3method(get_params,metab_bayes) S3method(get_params,metab_model) @@ -21,6 +22,7 @@ S3method(predict_DO,metab_Kmodel) S3method(predict_DO,metab_model) S3method(predict_DO,metab_night) S3method(predict_DO,metab_sim) +S3method(predict_metab,metab_2station) S3method(predict_metab,metab_Kmodel) S3method(predict_metab,metab_bayes) S3method(predict_metab,metab_model) @@ -97,6 +99,7 @@ export(revise) export(sim_Kb) export(sim_pred_Kb) export(specs) +exportClasses(metab_2station) exportClasses(metab_Kmodel) exportClasses(metab_bayes) exportClasses(metab_mle) diff --git a/R/data.R b/R/data.R new file mode 100644 index 0000000..cc45a05 --- /dev/null +++ b/R/data.R @@ -0,0 +1,44 @@ +#' Example two-station (VFTS) input data +#' +#' A 30-day example dataset for fitting a two-station (upstream/downstream, +#' a.k.a. VFTS) metabolism model with \code{\link{metab_2station}}. It is a +#' subset of the \code{VFTS-2} (variable-travel-time) run from the published +#' two-station metabolism modeling dataset for a reach of the Colorado River +#' in Glen Canyon, covering 2011-07-31 through 2011-08-29 at the source +#' data's native 15-minute timestep, plus a lead-in block of upstream DO +#' observations (exactly as long as the longest travel time in the dataset +#' requires -- see \code{\link{metab_2station}}'s "Two-station data +#' requirements" section) immediately before 2011-07-31. +#' +#' @format A data.frame with 2904 rows and the 9 columns expected by +#' \code{\link{metab_2station}}'s \code{data} argument, each carrying +#' \code{\link[unitted]{unitted}} units matching \code{\link{mm_data}}: +#' \describe{ +#' \item{solar.time}{POSIXct timestamp, UTC} +#' \item{DO.obs.up}{dissolved oxygen observed at the upstream station, +#' mg O2 / L} +#' \item{DO.sat.up}{dissolved oxygen at equilibrium saturation at the +#' upstream station, mg O2 / L} +#' \item{DO.obs.down}{dissolved oxygen observed at the downstream +#' station, mg O2 / L} +#' \item{DO.sat.down}{dissolved oxygen at equilibrium saturation at the +#' downstream station, mg O2 / L} +#' \item{light}{photosynthetically active radiation, umol m^-2 s^-1} +#' \item{depth}{reach depth, m} +#' \item{temp.water}{water temperature at the downstream station, degC} +#' \item{travel.time}{reach travel time between the upstream and +#' downstream stations, d} +#' } +#' +#' @source Filtered to \code{model_run == 'VFTS-2'}; see +#' \code{data-raw/two_station_example.R} for the extraction/renaming code. +#' +#' Bishop, I.W., Deemer, B.R., Kennedy, T.A., Payn, R.A., Hall Jr, R.O. and +#' Yackulic, C.B., 2026. A simplified two-station approach for modeling +#' metabolism in dam tailwaters subject to diel flow variation. Limnology +#' and Oceanography: Methods, p.e70066. +#' \url{https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lom3.70066} +#' +#' Data archived at ScienceBase: +#' \url{https://www.sciencebase.gov/catalog/item/6887d457d4be024722b4aae2} +"two_station_example" diff --git a/R/metab.R b/R/metab.R index 8a45c07..6e6b705 100644 --- a/R/metab.R +++ b/R/metab.R @@ -59,11 +59,12 @@ metab <- function(specs=specs(mm_name()), data=v(mm_data(NULL)), data_daily=v(mm model_type <- mm_parse_name(specs$model_name)$type metab_fun <- switch( model_type, - bayes = metab_bayes, - Kmodel = metab_Kmodel, - mle = metab_mle, - night = metab_night, - sim = metab_sim) + bayes = metab_bayes, + bayes_2station = metab_2station, + Kmodel = metab_Kmodel, + mle = metab_mle, + night = metab_night, + sim = metab_sim) # run the model metab_fun(specs=specs, data=data, data_daily=data_daily, info=info) diff --git a/R/metab_2station.R b/R/metab_2station.R index 36efca7..e298dd5 100644 --- a/R/metab_2station.R +++ b/R/metab_2station.R @@ -1,15 +1,26 @@ -#' Two-station Bayesian metabolism model fitting function (stub) +#' @include metab_model-class.R metab_bayes.R +NULL + +#' Two-station Bayesian metabolism model fitting function #' #' Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to -#' estimate GPP and ER from paired upstream and downstream DO, temperature, -#' light, and travel-time data. This function is currently a stub: it -#' validates the \code{data} argument and enforces two-station-specific data -#' requirements, but does not yet fit a model. See \code{\link{mm_name}} to -#' choose a Bayesian model and \code{\link{specs}} for relevant options for -#' the \code{specs} argument. +#' estimate GPP, ER, and K600 from paired upstream and downstream DO, +#' temperature, light, and travel-time data, using the single fixed Stan +#' model in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See +#' \code{\link{mm_name}} to choose a Bayesian model and \code{\link{specs}} +#' for relevant options for the \code{specs} argument. +#' +#' Unlike \code{\link{metab_bayes}}, which supports many model structures via +#' \code{split_dates}/\code{pool_K600}/etc., \code{metab_2station} always +#' fits every date jointly in a single Stan call (\code{specs$split_dates} is +#' forced to \code{FALSE} by \code{\link{specs}}), because the +#' upstream-downstream lag shift ties each date's first modeled rows to the +#' previous date's last rows. #' #' @inheritParams metab -#' @return Not yet implemented; currently always errors after data validation. +#' @return A metab_2station object containing the fitted model. This object +#' can be inspected with the functions in the +#' \code{\link{metab_model_interface}} and also \code{\link{get_mcmc}}. #' #' @section Two-station data requirements: In addition to the checks #' performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the @@ -30,57 +41,338 @@ metab_2station <- function( info=NULL ) { - # Check data for correct column names & units - dat_list <- mm_validate_data(data, data_daily, 'metab_2station') + stanfit <- NULL + fitting_time <- system.time({ + # Check data for correct column names & units + dat_list <- mm_validate_data(data, data_daily, 'metab_2station') + data_v <- v(dat_list$data) + + travel_time <- data_v$travel.time + solar_time <- data_v$solar.time + + # a. travel.time must be strictly positive + if(any(travel_time <= 0)) { + stop('travel.time must be > 0') + } + + # b. travel.time is expected in days, so any value >= 1 almost certainly + # reflects a units mistake (e.g., minutes or hours rather than days) + if(any(travel_time >= 1)) { + stop('travel.time must be < 1 day; values >= 1 suggest incorrect units (expected days)') + } + + # c. there must be enough lead-in rows of upstream DO before the first + # modeled row to cover the longest travel time in the dataset. timestep_days + # is the median observation interval, in days; max_lag is the number of + # timesteps by which upstream data must lead downstream predictions. The + # first max_lag rows of data serve only as lead-in and cannot themselves be + # modeled, so at least max_lag + 1 rows are required overall. + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + max_lag <- max(round(travel_time / timestep_days)) + if(nrow(dat_list$data) <= max_lag) { + lead_in_needed <- max_lag - nrow(dat_list$data) + 1 + stop(paste0( + 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', + max_lag, ' timestep(s), but only ', nrow(dat_list$data), ' row(s) were supplied; ', + 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) + } + + # Reconstruct the same "modeled rows" (post-lead-in-trim) index set that + # mm_ts_prep_data() computes internally, using the identical timestep_days/ + # lag/max_lag formula, so that Stan's date_index/time_index can be mapped + # back to actual dates and solar.times for the daily and instantaneous + # results below. (Duplicated here rather than exposed by + # mm_ts_prep_data(), which returns only the Stan-ready matrices.) + n_total <- nrow(dat_list$data) + keep <- seq.int(max_lag + 1, n_total) + modeled_solar_time <- solar_time[keep] + modeled_dates <- as.Date(modeled_solar_time) + date_df <- tibble::tibble(date=unique(modeled_dates), date_index=seq_along(unique(modeled_dates))) + n_days <- nrow(date_df) + if(length(keep) %% n_days != 0) { + stop(paste0( + 'dates have differing numbers of modeled rows after lead-in removal; ', + 'observations cannot be combined into a matrix: ', + paste(sprintf('%s (%d rows)', names(table(modeled_dates)), table(modeled_dates)), collapse=', '))) + } + n_obs <- length(keep) / n_days + obs_index_df <- tibble::tibble( + solar.time=modeled_solar_time, + DO.obs.down=data_v$DO.obs.down[keep], + date_index=rep(date_df$date_index, each=n_obs), + time_index=rep(seq_len(n_obs), times=n_days)) + + # Prepare the Stan data list (matrices from data, plus scalar priors from + # specs). modifyList (not c()) is used because mm_ts_prep_data() already + # supplies K600_lnorm_meanlog/K600_lnorm_sdlog (with its own fallback + # defaults), and those two names are also in specs$params_in; a plain + # c() would create duplicate-named list elements instead of overriding + data_list <- mm_ts_prep_data(dat_list$data, specs=specs) + data_list <- modifyList(data_list, specs[specs$params_in]) + + # Check and parse model file path + specs$model_path <- mm_locate_filename(specs$model_name) + + # determine how many cores to use, as in runstan_bayes() + tot_cores <- parallel::detectCores() + if(!is.finite(tot_cores)) tot_cores <- 1 + n_cores <- min(tot_cores, specs$n_cores) + + # Fit the model, collecting errors/warnings as strings rather than + # letting a bad dataset halt execution without reporting anything back + stop_strs <- character(0) + warn_strs <- character(0) + daily <- NULL + inst <- NULL + withCallingHandlers( + tryCatch({ + if(!suppressPackageStartupMessages(require(rstan))) { + stop("the rstan package is required for Stan MCMC models") + } + + consolelog <- utils::capture.output( + stanfit <- rstan::stan( + file=specs$model_path, data=data_list, pars=specs$params_out, + chains=specs$n_chains, cores=n_cores, + iter=specs$burnin_steps + specs$saved_steps, warmup=specs$burnin_steps, + thin=specs$thin_steps, verbose=specs$verbose, open_progress=FALSE), + split=specs$verbose) + + if(stanfit@mode == 2L) { + stop(paste(utils::capture.output(print(stanfit)), collapse='\n')) + } - travel_time <- v(dat_list$data$travel.time) - solar_time <- v(dat_list$data$solar.time) + # format the Stan summary matrix into per-variable data.frames + stan_mat <- rstan::summary(stanfit)$summary + mcmc_out <- format_mcmc_mat_nosplit( + stan_mat, data_list$n_days, data_list$n_obs, specs$model_name, + keep_mcmc=isTRUE(specs$keep_mcmcs), stanfit) - # a. travel.time must be strictly positive - if(any(travel_time <= 0)) { - stop('travel.time must be > 0') + # daily GPP/ER/K600 estimates: join Stan's date_index back to dates + date_index <- time_index <- index <- '.dplyr.var' + daily <- mcmc_out$daily %>% + dplyr::left_join(date_df, by='date_index') %>% + dplyr::select(-date_index, -time_index, -index) %>% + dplyr::select(date, dplyr::everything()) + + # instantaneous DO.mod.down estimates come from the 'metab' Stan + # transformed parameter (posterior median), which format_mcmc_mat_nosplit() + # buckets by row count rather than by name since 'metab' isn't in its + # par_homes lookup table; find that bucket by its column names instead + is_metab_bucket <- sapply(mcmc_out, function(df) is.data.frame(df) && any(grepl('^metab_', names(df)))) + metab_bucket_name <- names(mcmc_out)[is_metab_bucket][1] + if(is.na(metab_bucket_name)) { + stop("could not find 'metab' in the Stan output; check that specs$params_out includes 'metab'") + } + inst <- mcmc_out[[metab_bucket_name]] %>% + dplyr::select(date_index, time_index, DO.mod.down=metab_50pct) %>% + dplyr::inner_join(obs_index_df, by=c('date_index','time_index')) %>% + dplyr::select(solar.time, DO.obs.down, DO.mod.down) %>% + dplyr::arrange(solar.time) + + }, error=function(err) { + stop_strs <<- c(stop_strs, err$message) + }), warning=function(war) { + warn_strs <<- c(warn_strs, war$message) + invokeRestart("muffleWarning") + }) + + # if fitting failed, fill in NA daily estimates (with real dates) so the + # returned model at least reports which dates were attempted + if(length(stop_strs) > 0 || is.null(daily)) { + na_vec <- rep(as.numeric(NA), nrow(date_df)) + daily <- data.frame( + date=date_df$date, + GPP_daily_2.5pct=na_vec, GPP_daily_50pct=na_vec, GPP_daily_97.5pct=na_vec, + ER_daily_2.5pct=na_vec, ER_daily_50pct=na_vec, ER_daily_97.5pct=na_vec, + K600_daily_2.5pct=na_vec, K600_daily_50pct=na_vec, K600_daily_97.5pct=na_vec) + inst <- NULL + } + daily <- dplyr::mutate(daily, valid_day=TRUE, warnings='', errors='') + + fit <- list( + daily=daily, inst=inst, + warnings=trimws(unique(warn_strs)), errors=trimws(unique(stop_strs))) + }) + + # Package and return results + mm <- metab_model( + model_class="metab_2station", + info=info, + fit=fit, + log=NULL, + mcmc=if(isTRUE(specs$keep_mcmcs)) stanfit else NULL, + mcmc_data=if(isTRUE(specs$keep_mcmc_data)) data_list else NULL, + fitting_time=fitting_time, + compile_time=system.time({}), # rstan::stan() compiles & samples in one call; not timed separately + specs=specs, + data=dat_list$data, # keep the units if given + data_daily=dat_list$data_daily) + + # Update data with DO predictions + success <- !is.null(fit$inst) && length(fit$errors) == 0 + if(success) { + mm@data <- predict_DO(mm) + } else { + warntxt <- paste0( + 'Modeling failed\n', + if(length(fit$warnings) > 0) paste0(' Warnings:\n', paste0(' ', fit$warnings, collapse='\n')), + if(length(fit$errors) > 0) paste0(' Errors:\n', paste0(' ', fit$errors, collapse='\n'))) + warning(warntxt) } - # b. travel.time is expected in days, so any value >= 1 almost certainly - # reflects a units mistake (e.g., minutes or hours rather than days) - if(any(travel_time >= 1)) { - stop('travel.time must be < 1 day; values >= 1 suggest incorrect units (expected days)') + # Return + mm +} + + +#### metab_2station class #### + +#' Metabolism model fitted by two-station (VFTS) Bayesian MCMC +#' +#' \code{metab_2station} models use Bayesian MCMC methods to fit values of +#' GPP, ER, and K600 from paired upstream/downstream DO curves. This class +#' inherits from \code{metab_bayes} (same \code{log}/\code{mcmc}/ +#' \code{mcmc_data}/\code{compile_time} slots, and therefore the same +#' \code{\link{get_mcmc}}, \code{\link{get_mcmc_data}}, and +#' \code{\link{get_log}} methods), but \code{predict_metab} and +#' \code{predict_DO} are overridden below because two-station's fitted-value +#' structure and output columns differ from one-station's. +#' +#' @exportClass metab_2station +#' @family metab.model.classes +setClass("metab_2station", contains="metab_bayes") + + +#' @describeIn get_params Does the same Stan-output-to-streamMetabolizer +#' renaming as \code{get_params.metab_bayes}, but (unlike that method) does +#' not delegate the rest of the work to \code{get_params.metab_model} via +#' \code{NextMethod()}: that generic implementation looks up parameter +#' names via \code{get_param_names()}, which assumes a +#' \code{metab_()}-named model constructor (would look for +#' \code{metab_bayes_2station}, which doesn't exist -- our constructor is +#' \code{metab_2station}) and streamMetabolizer's ODE-based dDOdt framework +#' for one-station models, neither of which apply to the two-station +#' steady-state model. \code{fixed} column/star annotations (relevant only +#' to models that can take fixed daily parameters from \code{data_daily}) +#' are not supported here. +#' @export +#' @import dplyr +get_params.metab_2station <- function( + metab_model, date_start=NA, date_end=NA, uncertainty=c('sd','ci','none'), messages=TRUE, ...) { + + uncertainty <- match.arg(uncertainty) + + fit <- metab_model@fit$daily + if(is.null(fit)) return(NULL) + + # Stan prohibits '.' in variable names, so convert back from '_' to '.', + # as in get_params.metab_bayes + parnames <- setNames(gsub('_', '\\.', metab_model@specs$params_out), metab_model@specs$params_out) + parnames <- parnames[order(nchar(parnames), decreasing=TRUE)] + for(i in seq_along(parnames)) { + names(fit) <- gsub(names(parnames[i]), parnames[[i]], names(fit)) } + names(fit) <- gsub('_mean$', '', names(fit)) + names(fit) <- gsub('_sd$', '.sd', names(fit)) + names(fit) <- gsub('_50pct$', '.median', names(fit)) + names(fit) <- gsub('_2.5pct$', '.lower', names(fit)) + names(fit) <- gsub('_97.5pct$', '.upper', names(fit)) + + fit <- mm_filter_dates(fit, date_start=date_start, date_end=date_end) - # c. there must be enough lead-in rows of upstream DO before the first - # modeled row to cover the longest travel time in the dataset. timestep_days - # is the median observation interval, in days; max_lag is the number of - # timesteps by which upstream data must lead downstream predictions. The - # first max_lag rows of data serve only as lead-in and cannot themselves be - # modeled, so at least max_lag + 1 rows are required overall. - timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) - max_lag <- max(round(travel_time / timestep_days)) - if(nrow(dat_list$data) <= max_lag) { - lead_in_needed <- max_lag - nrow(dat_list$data) + 1 - stop(paste0( - 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', - max_lag, ' timestep(s), but only ', nrow(dat_list$data), ' row(s) were supplied; ', - 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) + metab.vars <- c('GPP.daily', 'ER.daily', 'K600.daily') + for(mv in metab.vars) { + if(paste0(mv, '.median') %in% names(fit)) fit[[mv]] <- fit[[paste0(mv, '.median')]] } + keep.cols <- c('date', unlist(lapply(metab.vars, function(mv) grep(paste0('^', mv, '($|\\.)'), names(fit), value=TRUE)))) + params <- fit[intersect(keep.cols, names(fit))] + + params <- switch( + uncertainty, + 'none' = params[!grepl('\\.median$|\\.sd$|\\.lower$|\\.upper$', names(params))], + 'sd' = params[!grepl('\\.median$|\\.lower$|\\.upper$', names(params))], + 'ci' = params[!grepl('\\.median$|\\.sd$', names(params))]) - stop('metab_2station not yet implemented') + # attach raw warnings/errors columns (not yet compressed into a single + # column); show()'s pretty_print_ddat()/compress_msgs() does that + # compression itself at print time, as in get_params.metab_model + if(messages && exists('date', fit) && any(c('warnings','errors') %in% names(fit))) { + msgs <- fit[c('date','warnings','errors') %>% { .[. %in% names(fit)] }] + params <- left_join(params, msgs, by='date', copy=TRUE) + } + + params } -#' @describeIn predict_DO Stub for two-station (VFTS) models. Not yet -#' functional: \code{metab_2station} does not yet produce a fitted Stan -#' model to predict from (see \code{\link{metab_2station}}), so this always -#' errors. Once implemented, this method will return a data.frame + +#' @describeIn predict_metab Pulls daily GPP, ER, and K600 estimates out of +#' the two-station Stan model results. +#' @export +#' @import dplyr +predict_metab.metab_2station <- function(metab_model, date_start=NA, date_end=NA, ...) { + + Var1 <- Var2 <- '.dplyr.var' + fit.names <- expand.grid(c('50pct','2.5pct','97.5pct'), c('GPP_daily','ER_daily','K600_daily'), stringsAsFactors=FALSE) %>% + select(Var2, Var1) %>% + apply(MARGIN=1, FUN=function(row) do.call(paste, c(as.list(row), list(sep='_')))) + metab.names <- expand.grid(c('','.lower','.upper'), c('GPP','ER','K600'), stringsAsFactors=FALSE) %>% + select(Var2, Var1) %>% + apply(MARGIN=1, FUN=function(row) do.call(paste0, as.list(row))) + + fit <- metab_model@fit$daily %>% + mm_filter_dates(date_start=date_start, date_end=date_end) + if(is.null(fit) || !all(fit.names %in% names(fit))) { + stop('could not find GPP_daily, ER_daily, and K600_daily estimates in the model fit') + } + preds <- fit[c('date', fit.names)] %>% + setNames(c('date', metab.names)) + + # add date-specific fitting warnings/errors, as in predict_metab.metab_bayes + warnings <- errors <- '.dplyr.var' + if(!is.null(fit) && all(c('date','warnings','errors') %in% names(fit))) { + messages <- fit %>% + select(date, warnings, errors) %>% + compress_msgs('msgs.fit', warnings.overall=metab_model@fit$warnings, errors.overall=metab_model@fit$errors) + preds <- full_join(preds, messages, by='date', copy=TRUE) + } else { + preds <- mutate(preds, msgs.fit=NA) + } + + preds <- mutate( + preds, + warnings=if(length(metab_model@fit$errors) > 0) NA else '', + errors=if(length(metab_model@fit$errors) > 0) NA else '') + + preds +} + + +#' @describeIn predict_DO Two-station (VFTS) models. Returns a data.frame #' with columns \code{solar.time}, \code{DO.obs.down} (the observed #' downstream DO from the input data), and \code{DO.mod.down} (the -#' two-station Stan model's predicted downstream DO) -- unlike the -#' one-station \code{predict_DO} methods, which return \code{DO.obs}/ -#' \code{DO.mod}. +#' posterior median of the two-station Stan model's fitted downstream DO) +#' -- unlike the one-station \code{predict_DO} methods, which return +#' \code{DO.obs}/\code{DO.mod}. The values are those computed once at +#' fitting time (see \code{\link{metab_2station}}); \code{use_saved=FALSE} +#' (on-demand recomputation from the fitted daily GPP/ER/K600 medians) is +#' not implemented. #' @export predict_DO.metab_2station <- function(metab_model, date_start=NA, date_end=NA, ..., use_saved=TRUE) { + if(!isTRUE(use_saved)) { + stop("predict_DO(use_saved=FALSE) is not implemented for metab_2station; only the fitted-time DO.mod.down values are available") + } + + inst <- metab_model@fit$inst + if(is.null(inst)) { + stop("no DO.mod.down predictions are available; the model fit may have failed (see get_fit(metab_model))") + } + # NOTE: R/plot_DO_preds.R and tests/testthat/helper-rmse_DO.R both - # hard-code the one-station DO.obs/DO.mod column names; they'll need to - # branch on (or be parameterized for) DO.obs.down/DO.mod.down once this - # method actually returns two-station predictions. - stop("predict_DO for two-station models not yet implemented — requires completion") + # hard-code the one-station DO.obs/DO.mod column names; they still need to + # branch on (or be parameterized for) DO.obs.down/DO.mod.down before those + # tools will work with two-station predictions -- deferred, out of scope + # for this PR. + mm_filter_dates(inst, date_start=date_start, date_end=date_end) } diff --git a/R/mm_name.R b/R/mm_name.R index cc41d30..4df02fc 100644 --- a/R/mm_name.R +++ b/R/mm_name.R @@ -54,6 +54,9 @@ #' @param type character. The model type. Options: \itemize{ \item \code{mle}: #' maximum likelihood estimation (see also \code{\link{metab_mle}}) \item #' \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item +#' \code{bayes_2station}: two-station (upstream/downstream, a.k.a. VFTS) +#' Bayesian model with a single fixed structure (see also +#' \code{\link{metab_2station}}) \item #' \code{night}: nighttime regression (see also \code{\link{metab_night}}) #' \item \code{Kmodel}: regression of \emph{daily} estimates of #' \code{K600.daily} versus discharge, time, etc., usually for 3-phase @@ -155,7 +158,7 @@ #' mm_name('sim', err_proc_acor=TRUE) #' mm_name('bayes', pool_K600='binned') mm_name <- function( - type=c('mle','bayes','night','Kmodel','sim'), + type=c('mle','bayes','bayes_2station','night','Kmodel','sim'), #pool_GPP='none', pool_ER='none', pool_eoi='alldays', pool_epc='alldays', pool_epi='alldays', pool_K600=c('none', 'normal','normal_sdzero','normal_sdfixed', @@ -175,9 +178,31 @@ mm_name <- function( engine=c('stan','nlm','lm','mean','loess','rnorm'), check_validity=TRUE) { - # determine type - type <- match.arg(type) - + # determine type. 'bayes_2station' is matched exactly, before match.arg's + # partial-prefix matching, because 'b' would otherwise be an ambiguous + # abbreviation between 'bayes' and 'bayes_2station' -- so unlike the other + # types, 'bayes_2station' must be spelled out in full (no abbreviations). + # match.arg's choices are narrowed to exclude it so that pre-existing + # abbreviations like 'b' (-> 'bayes') and 'm' (-> 'mle') stay unambiguous. + if(missing(type)) { + type <- eval(formals(mm_name)$type)[1] + } else if(length(type) == 1 && identical(type, 'bayes_2station')) { + type <- 'bayes_2station' + } else { + type <- match.arg(type, choices=setdiff(eval(formals(mm_name)$type), 'bayes_2station')) + } + + # bayes_2station has a single fixed model structure rather than being built + # from combinations of pool_K600/err_*/ode_method/GPP_fun/ER_fun/ + # deficit_src/engine, so skip the argument-combination machinery below and + # return the one valid name directly + if(type == 'bayes_2station') { + mmname <- 'b2_np_oi_tr_plrckm.stan' + check_validity <- if(!is.logical(check_validity)) stop("need check_validity to be a logical of length 1") else check_validity[1] + if(isTRUE(check_validity)) mm_validate_name(mmname) + return(mmname) + } + # set type-specific defaults where values weren't specified . <- '.dplyr.var' if(type != 'Kmodel') { diff --git a/R/mm_valid_names.R b/R/mm_valid_names.R index 69a7111..1a98d1b 100644 --- a/R/mm_valid_names.R +++ b/R/mm_valid_names.R @@ -10,7 +10,7 @@ #' @examples #' mm_valid_names('mle') #' @export -mm_valid_names <- function(type=c('bayes','mle','night','Kmodel','sim')) { +mm_valid_names <- function(type=c('bayes','bayes_2station','mle','night','Kmodel','sim')) { type <- match.arg(type, several.ok=TRUE) @@ -39,6 +39,11 @@ mm_valid_names <- function(type=c('bayes','mle','night','Kmodel','sim')) { mnames <- grep('^b_', dir(system.file('models', package='streamMetabolizer')), value=TRUE) favorites <- c('b_np_oipi_tr_plrckm.stan','b_np_oi_tr_plrckm.stan','b_np_pi_tr_plrckm.stan','b_np_oipp_tr_plrckm.stan') }, + bayes_2station={ + # single fixed model structure; no combinatorial name-building needed + mnames <- 'b2_np_oi_tr_plrckm.stan' + favorites <- mnames + }, mle={ opts <- expand.grid( type='mle', diff --git a/R/specs.R b/R/specs.R index 637c1a5..a1facb8 100644 --- a/R/specs.R +++ b/R/specs.R @@ -250,6 +250,12 @@ #' proportional to light (with noise) and is applied to GPP rather than to #' dDO/dt. #' +#' @param K600_lnorm_meanlog hyperparameter for \code{type='bayes_2station'}. +#' The mean of a lognormal prior distribution for K600_daily. +#' @param K600_lnorm_sdlog hyperparameter for \code{type='bayes_2station'}. The +#' standard deviation parameter of a lognormal prior distribution for +#' K600_daily. +#' #' @param params_in Character vector of hyperparameters to pass from the specs #' list into the data list for the MCMC run. Will be automatically generated #' during the specs() call; need only be revised if you're using a custom @@ -423,7 +429,12 @@ specs <- function( err_proc_acor_phi_beta = 1, err_proc_acor_sigma_scale = 1, err_mult_GPP_sdlog_sigma = 1, - + + # hyperparameters for two-station (bayes_2station) K600. GPP_daily_mu, + # GPP_daily_sigma, ER_daily_mu, and ER_daily_sigma above are reused as-is + K600_lnorm_meanlog = 2.484907, + K600_lnorm_sdlog = 1.0, + # vector of hyperparameters to include as MCMC data params_in, @@ -624,9 +635,63 @@ specs <- function( warning(e) return(model_name) }) - if(features$engine == "NA") + if(features$engine == "NA") stop('engine must be specified for Bayesian models') - + + }, + 'bayes_2station' = { + + # bayes_2station has a single fixed model structure (see + # inst/models/b2_np_oi_tr_plrckm.stan), so params_in/params_out are + # hardcoded here rather than built up from pool_K600/GPP_fun/ER_fun/ + # err_* toggles as in the 'bayes' case above + + # the six scalar prior hyperparameters spliced into the Stan data list + # by mm_ts_prep_data() + all_specs$params_in <- c( + 'GPP_daily_mu', 'GPP_daily_sigma', + 'ER_daily_mu', 'ER_daily_sigma', + 'K600_lnorm_meanlog', 'K600_lnorm_sdlog') + + # list all needed arguments + included <- c( + # model setup + 'model_name', 'engine', 'split_dates', 'keep_mcmcs', 'keep_mcmc_data', + + # params_in is both a vector of specs to include and a vector to include in specs + all_specs$params_in, 'params_in', + + # inheritParams runstan_bayes + 'params_out', 'n_chains', 'n_cores', + 'burnin_steps', 'saved_steps', 'thin_steps', 'verbose' + ) + + # compute some arguments + if('engine' %in% yes_missing) { + all_specs$engine <- 'stan' + } + if('split_dates' %in% yes_missing) { + # forced FALSE: the upstream/downstream lag shift ties consecutive + # days together, so days can't be modeled independently + all_specs$split_dates <- FALSE + } + if('params_out' %in% yes_missing) { + # 'metab' (the Stan transformed-parameter matrix of modeled + # downstream DO) is included so that predict_DO() has a fitted value + # to report; without it the model would only yield daily GPP/ER/K600 + all_specs$params_out <- c('GPP_daily', 'ER_daily', 'K600_daily', 'sigma', 'metab') + } + + # check for errors/inconsistencies + model_path <- tryCatch( + mm_locate_filename(model_name), + error=function(e) { + warning(e) + return(model_name) + }) + if(features$engine == "NA") + stop('engine must be specified for Bayesian models') + }, 'mle' = { # determine which init values will be needed diff --git a/data-raw/two_station_example.R b/data-raw/two_station_example.R new file mode 100644 index 0000000..0142b89 --- /dev/null +++ b/data-raw/two_station_example.R @@ -0,0 +1,89 @@ +# Builds data/two_station_example.rda from the VFTS paper's published input +# data. Not run automatically as part of the package build/check; run +# manually (with the package root as the working directory) whenever the +# example dataset needs to be regenerated. +# +# File source: 2_station/Data/2_VFTS_and_One-station_model_input.csv, a sibling +# directory of the package root (not included in the package itself) + +# Download from ScienceBase: https://www.sciencebase.gov/catalog/item/6887d457d4be024722b4aae2 + +# Paper URL: https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lom3.70066 + +library(dplyr) +library(unitted) + +# --- read & filter ----------------------------------------------------- + +raw <- read.csv( + file.path('..', '2_station', 'Data', '2_VFTS_and_One-station_model_input.csv'), + stringsAsFactors=FALSE) + +# VFTS-2 is the variable-travel-time two-station run (as opposed to the +# fixed-travel-time VFTS-1/VFTS-3 runs or the one-station-only 'OS' run) and +# covers an intermediate date range (2008-03-11 to 2014-02-27) +vfts2 <- raw %>% + filter(model_run == 'VFTS-2') %>% + mutate(datetime = as.POSIXct(datetime, format='%Y-%m-%dT%H:%M:%SZ', tz='UTC')) %>% + arrange(datetime) + +# 30 consecutive days from the middle of the dataset (2011, roughly halfway +# between 2008 and 2014), avoiding the start/end edges of both the overall +# dataset and of this particular gap-free stretch of observations (verified +# separately to run gap-free from 2011-07-28 to 2011-08-30 at the source +# data's native 15-minute timestep). +modeled_start <- as.POSIXct('2011-07-31 00:00:00', tz='UTC') +modeled_end <- as.POSIXct('2011-08-29 23:45:00', tz='UTC') +timestep_days <- 15/(24*60) + +# metab_2station()'s upstream-DO lag shift (see mm_ts_prep_data() and +# metab_2station()'s "Two-station data requirements" section) needs +# max_lag = max(round(travel.time / timestep_days)) rows of lead-in +# immediately before modeled_start -- and because mm_ts_prep_data() trims +# max_lag rows off the *start of the whole array*, not off each calendar +# day, that lead-in window must be exactly max_lag rows (not e.g. a whole +# extra day) or the first modeled date ends up with a different row count +# than the rest, which mm_ts_prep_data() rejects. max_lag is computed from a +# generous 2-day candidate lead-in window and then trimmed to size. +candidate_start <- modeled_start - as.difftime(2, units='days') +candidate <- vfts2 %>% filter(datetime >= candidate_start, datetime <= modeled_end) +max_lag <- max(round(candidate$travel_time / timestep_days)) +lead_in_start <- modeled_start - as.difftime(max_lag * 15, units='mins') + +vfts2_window <- vfts2 %>% + filter(datetime >= lead_in_start, datetime <= modeled_end) + +# confirm the window is gap-free at the native 15-min timestep, and that +# trimming the lead-in rows (as mm_ts_prep_data() does) leaves exactly 30 +# modeled dates with equal row counts +stopifnot(all(abs(diff(as.numeric(vfts2_window$datetime)) - 15*60) < 1e-6)) +modeled_dates <- as.Date(vfts2_window$datetime[seq.int(max_lag+1, nrow(vfts2_window))]) +stopifnot(length(unique(table(modeled_dates))) == 1) +stopifnot(length(unique(modeled_dates)) == 30) + +# --- rename to package conventions & attach units ----------------------- + +renamed <- vfts2_window %>% + transmute( + solar.time = datetime, + DO.obs.up = upstream_DO, + DO.sat.up = upstream_DO_sat, + DO.obs.down = downstream_DO, + DO.sat.down = downstream_DO_sat, + light = light, + depth = reach_depth, + temp.water = downstream_temp, + travel.time = travel_time) +# model_run and lag are intentionally dropped (not package columns) + +template <- mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, + light, depth, temp.water, travel.time) +two_station_example <- renamed +for(col in names(template)) { + two_station_example[[col]] <- u(renamed[[col]], get_units(template[[col]])) +} +two_station_example <- two_station_example[names(template)] + +# --- save ----------------------------------------------------------------- + +usethis::use_data(two_station_example, overwrite=TRUE) diff --git a/data/two_station_example.rda b/data/two_station_example.rda new file mode 100644 index 0000000000000000000000000000000000000000..d6292ef0795293d2ffb132683667981a42385aa6 GIT binary patch literal 17629 zcmeIZWn3Ih@Ggoh2@rxyAV6So4;nlW77Y+IxP;)a5M)^_1W90Vhv06DyK8U_9^Bns z637|y{?9qT`{{nTUvII@?DWh|cXdBqRXy|68(m8wZb=fxoCb#Y38}(kPonW?{mFAm*uY8*R$1m+_)|R~=&fue<{(RV@rlYEO)A%rB zlHg0C`A7gK(jSWd398(Gd-TW6#xp+n1bZ)JE*iFv7(r5rLD)}T;eOI@oAERAvq7y> z-2T;;*lyi!^N3sl#gYQ=Ld8pc5(UK)`29mImXEx7EJyS5?^IAw+(UdKGj1dB@pqFH*6!i%`Lj}Y z8l3`(;ge5aPrkcQANLp^=ON)|$-Yy)ckM*d3rx7D{Y4`E6cH$2F<5=FWD-!*|I-=; z<&pP)-~Mk2ih@^;JT@fn*55xhmckMz1pg8Z^@8oyUB~V^#JkCdX3Jlll$=7`vk0dVe24krcy)FN0 zqQw*t_vanS2SLSl-lP@KzwY;nlury;&n2SqV13-AJ|SBWBgPcvm6u5~ekU*XZjy(` zo^O(=;5OnPOzMFEr-vhzk(t6MtE%=PvHH z&<8d)pu!R;N`{pubW=h5HUf^k<``Xetff_~EkA#}f?bxfHJqM}h``pGw30hSr@R%2 z8wM2Q8tWPgh4N0QO)DT@l-0JDX z&of9-En_|5^e&ZXYx@@a=w0&tea6+*Tb{HhYpm-pzG{bj(5d}dhIQuG7mli@mSir^ z&`d7STvOa4QzG&RS}lk2`2w`=_%Yt?CKX%Y*!&TVe!Sq}i-JWv ze%-Xa$v;aLty?tr6tLaZ|NpoDpL;;?`v1ti+57X?zf;vU%)}olPq?DQSIdCeexq|F zNkvyf!(@BrCO5wXjI4j=>RP!;eaCF)STD8f=HYgcSgiAc->vvq zIa_~7tmeDCel2?2KmLAp2xVoK%e;l=+JVVsnGb)c^`RE?yVd54L^EDz2=KEuBk))V)yn@E7tYDT2Gmm$X^$_Kc+c0KUFC~`JK&HrbH%gB$uX=iV&h$OM z$-Yi2wMMwdGOhX}`7|6BQbfY)8-tT&owWzceUQ~3c}WX?&I&(s7W#AlO;Y*G>kMK} z3*MRTbw>}e)Deq0mxo>Lv0qwHsnqv?%I)0?Z->!m9h^~wtI{+@dcL*9Ci9WS0@ex%%1i+|BtRmy}xQP<2f=Bke-X%7j44{$WL};IQpfz!S*ab7;*$(dLwm zt?5*CRdplJ(D2+x+ixZ-iK9|*}$`%d6<72Q0$TD-t(p zRkUB3Jy%g#F+6U$Il!F~Xk)6wETB_qmM2gZS_S*Y!Ec30OqG(1D&vA~1I>G4TXoR? zay+CT*UD}WJ#M7U;Ko|*E4jcI+gNpw;TuN^R8^$P}#C5?z;y$YAT ziUu~0(V2<$!`xj?hyj63txM#G9|0=fKLX^wR_yI~3(p$Z)#?>b?ZF!6uKH>lKp4w! zPEaWhtP-tqlicl2YX!Q8QgMe;?ECg-DlJ&Q9g>$mxx9YW?utNUM=1-W(=x0Q%F{;S@#sM zv$Nachq;hA$7*VtOoB0f#h3rucj3CL{kV83f?Zfum<_?W(A2KWrty19>FT-`CY~v{ z5}@h->;v#=x7O5Dk;Kqw(yw>i7>MU0;!^j0_FyMeul)6F>UMQ zPGH0C`(9+O@Qw7UZE+)lc3|W5=JS(^~Lv-ty#p#$Io$gn3 zLZothcMaI>A<#i+w&6MFqFEZgqCd}&=ey83;O z2kf`c;c⪻eD-5_3YVwII#V7-l+Z3c-6Sw??`dx2O!N48?rs}goK2$nK2M&!(Peu z0{rf<7H`l-FBb(2s@}<+0DDw5N}VDaRZ>l;km3AP1;T@U z9;#DBy=1x@$oNI(&iCqsy=6s%MLeB|hR@<@?=l7fxFp&g!5bG2{|ORVNd8;HPF0#% z^`fPA>dvTWs=CR%y4aqe+>5vV&?#m2w)W1bL32V=(-4g8BE$Df-Vz!fiRaYRG)PB7 z?&SICn$%p{0+2gxJ()j*ixXQ7y#LPeGWYvao~JLPQ+%s-#|vR5Ff*7N0rd#ISSmPq z9&4yiwtlkGd?-o$cm3tvuv(rhb_%6NgOK!ShjJK$0cGqf2h>coK!}FAI!KljCA8+0 zbk5#rdZ{&kVy(ZFh_JZGw@ES5RfoYf0d&GmoXW^5s`%!FTVuZcc>N?FIO{7O2_3VvxB>S28(Rfu(< z?g^WrhLUTvA@OvKUII&f5T=d=3vkD-_WtGeuyjoV>~xe#66?UHZbXJ^HvJr5Zeoel zjAXJ_ZU(BN_l#Z2gT$~zG+g41u%ckYZk)}O{j7gP7=k)8f(bi$$<>*=b(sB3E`|RaMG!khPBsfD6M5O3l)gMIhL5{Cqfx=jguaiy@$uU?f2#7e`@3Qst!r_(k zxDwx}rfM=-9Z>kgZ*n9}2*6mv@<;Bqsh;o`hxZ3@+UD@V!I5q-?Ktf?NW!@n;f}#L zEb6-26TUZN76(`jdHN(qTa)8YpFY`s@>+<~F|p1BgTXw|zT8C1W*c^3ErF63p>QrGzH74N}6#nJWq7mMq24FCu{fT_JESDEK>Nc#9^Khk2aU|BzL$ zuE1_XzA;edP*Z6>&~8JuF;ZzzZocnHt=&Q+yxL-GYio@LgIkjY;iD*Rg~S{S!i?y6XV*%1$oPVvF@< zYCiI~aPe5!7c`3)9K4v3CU&&UMA$ld~~y# zf39CyI6V0yLBKakKtK@kU@=@*aIZX?r)eb6P+Srh_`ZsZ(&aFj$o^8>^LPlc7|YFl2i>(uPoTDD@!%)E*My7F-iBpWHoCJ!`{Ty zDi3dk&X@`fE5!C_5d^aKscENPF6WF0tGlE-h_cHtB$6ykW88$1$$nrJUqN}#rY!KkPNC;0{59q2ayc@k5dD2?{d? zXx)TtzWqkp;*Is@V8_}lxl5D?CwsW%5ewtxh2btv3l~x!RA@=Y9y`)bYU%?do2@iFC;3ct z`fG`!Q?f$XUXS8}i7Zub2G)EEGXzfC=SU#09)Q{nfY(=0lF6)Om7(UVo*nd?z7L7bS`E9{4zf>G>+bm{0>>DxQi^*yxcuN=+MScoNJP?aGL zEMP}xF#5MXGJG{jLTO65Rh6xvdHKGdG5R`@3i>Bw@}w?IE@tM=kG<3l-Qp=>fjzJk zmz0zbA9P4#g+`@Qre=dGvC#&E2P-Q?G8;K7hO?nAW&?KnnZp(801(y!t97KwLn{vD zGNllCl(^N!o^a0~N0?`L+~Sl6jo@Rd(Zj>}1}3xQs5m=FU3zpBEsLOUF!#WEY$Rp) zvut5mp$@aaAHIQo9f6(tYng~ooZu8124XB`=Ck*X`na1G&Svpaks>yF7D?5zT8Rs0 z#f9{#haF9@ou*EVIOmjw%@F6=pkZEN^kPo%IHHj~J0fl>q7D*cr58l6E>fE$FDuy* zM&#{_EnF$m6)!IYGUUOH&eaW)%3wZv9z#XT6m}qt2unZuVrWRSm4#-B@gts=79`i0 zsj4n5^y%BDm0liD6d@Kdp`_|D0ff*vN@sKR%j@V-AwnM>1m1hAavao`Y@uV!27Qny zq&Mh!CTJkQRFw|~X9TOFl&eHqKQhmX*ujs2`szv3fA?&&`;gGu#scJl`oZw0s{V(l z#*Z+cxtj{HQ^xslg%*7aH2@(B^`lVfD+$|)t_uB+35hws*cJ0WWmFV<#Pl)K+tl(A zX8DJGZ?(tTiKyL-vlQGkd|2wV1MFrM1T0`w31&R7Lfm_!-t2YA*T{B7mN$vrW)viSYt81MDCF+()fvQ*bPD$@vXo5VxyW!BzIBfN6>E3Gcz18 ztHeI}@om%?hPjA4($coi37*m(w{mzPjhY?_5<~=clK9nOg~SzN$#2MYyFioG)0h3=3f&;Tan3@tW~QA%HVD1AYBvuV>__8k#bie3YE< zbpRH^h5apFl*k{&>}9Zz8E(P~xAJxwJ?59dNZ-ioUa3_?ByJQa`w&8NoG7kDXYolA z|2qSdBYHMwq;Vivn%ZE2P(?mpTj-HB6-YbwL?b=IN>DB;>b>)N)`_wg2Ua_T9E|Ts z$T@TpjlCI`KtB31pt7Et0lZ4sq>vqrp_ZrpjA08azLy!L!Uy(l#bwKo0*9tc{ z7Iz)rTuGl0O#J4NV-zS7qhA^v1tQe1OhA9F#<{VnO&=_h!LO&RA|a8UqfX4d(*dun z_YRVskYD!I?&BjRLEwHTRxcH=IeXMc;rrBfl~z=0Rn66b3GahCL-ix`HBTmUFd_Ps}pu61xWrH<(>B_yiFD(6lBfyVm&n{Wbfk-u5O2@ zB0C3Tk_9=erqm;Kz?5JtVp+nZRL+9MACsnra=wIQc(F`h-+oggnKVU`U>WOlpu|a2 zl4dxE*;0SyzyhJ&VZ|4&NR)_=#imwU-zAUf)kM>;8U%K9vw3y&^UsReQp(=2O$owt& z=v*$upv17`QgNNPUV4l*yKBPm3Mo+Un36Y{Wb?I-`c#Rar8l!@x?4Udnbe9+zwoPcv-TYe=Fl^VqxNU z5=ux1fqZ>w;*8}^((A_d>*5kSg>d9RjEaugA2Fr+X=;R;5q#w`U*q5`DWHyQCP=-~ zldhw;?o^jEtJd2@zUQBggHX|M^$Br|=seMd!gPPQ+eEWjJ~rOKmj}h9P~wzaRaqHJ|;8Kp4ns{su}pPr7Oy=G@b1lPZTXi0c4c{KkQCV zVuz4dYr_BcP;bYh1cW6$GeQWKym|qVQ5C9Z*gze?x6R(c_cLy!QBKs_%{O4fBDBS- z^!DKlCZipRwL#P97P=9g?^-2HWE_gDLLN$ zPCFhCqHl{*dqWGML121kVILoJv-c$44mZhfjIQ~MNcPe$R_?67RigI=ooWlxle#m( zs#K{MM&wQ#;>Gk297YfDS-Pn0+uiV=+ThT%VfM`n&ILF=YbQrt57}f_eYv}ODv!y@ zobGy5su9zLO~=`e(PM$%g_sD(Z)Dk!YJ)c`b6{drXs1Fl0{giL6HWpfYt_bc3aMtuxu;?$;d zIuYn+W*=Cxclt?J!E~8hg)Hs81i3KPPij%3B3XR* zJ0*`PU!%QQNBnYJithV7sSLemdDqWoX-5qC5|T4ra~r{WLok_>k<*Z|!g@{AAxZ zAdrM)gd|a4nm&EiZLsIxn5z~?BN1)?+l8&~oE6C;%6mK&wXI1Vv$hpHw=N``3v-sy z(HR2M&&_AAo0UFG?`>PuxL=-)*gL=*Yq{VIQz^lAW4I8=s)VC=|NNef;`efsUP4Uz z2zeCWC-^$0YxWBt18+~z;&dILrySCKo6Vnu^lNEiX#ZTE!u%AWPLS&Moj*0d3n8yw zrp0^g!LwROHiQCG8JN0NJLqt;TKOCuRCsjEKCl=yUUQ|N6_vti>eNUam1HVKu_mX5 zxrg0G(4TlHXiGck9^ed}$}DVcG5v6UrR4khqcEhL1>5CFas75|!nvlHinDr}>}5;d z)`3mQjV9G+60-!7t~M3m!KA}93FRlaD%@ed^D8C;p#%g^xv@LOsV~sq7Yn2Wk_bjb zjxi(WXk_ujs1{;%C{>bU$zUlnKV_-9oK*>8@ILNpT$n# zP1lxAaZ_9Ll~j`Qs7htbZz7oUy_p2V!j_fi1r%M@qJ5+z*yU-8T)sW7Xm{jrbZs%# zcJy?1^5>Xxuknw1oCSp?6Op^0gvN6oMTlJ3elS>|xHLf(Rr`^6EAIb%*n(h;z;{Iesc;aUo@`PKObmtpyN z!+IvyvvTSYyQ%WEppyb({%&X>qr}(mzN*Y%IZ67qIEFyEk~wE*RjP*B$*P+7QwPb* zQIez*6|{uIHb;(2*x_18xd=ik}~}fl`q!D+dBXS%}V%O~{^!ePN5= z^v~PYtb+p&&#l(J$x}hMDricQh>*aXl>^P!NE}2j%>i0f;8C(3d3*q_4`EM;y|fFR z43+uxt9>(TT1XsjbtUKGhfaiA-e+2kO;smfqgp5Ce1 zsmIKS%%qkSl1{01E=Fmvk8}nmtSaYoAt{?6yx4leoN=#ws!3N&IVE!Bo;@9Wr*UduE`H*_461k}cVMjA}o6NGt%i*WD z@X>g#SbneBC(*S}w7*b0Z3T$!6g(fV{&AUx{PvMLZ>~O7Uk~fgnu+PR0G+`uthO5K z&N4I!a0L6ah_)EXIlDO}J2q$^jo zdF?Tw$}h@v+NYCnc}mQ~4eebmC}(32X4T@=c|8x~I}ssr)2;^1MC~&t&re3EIe+iv z>m7x`6BXVM<2oA$xfBufq}HaS$ECS-KY{kI6Ut>2*Tdn{JUL3LxSTG<#d3LFy`m0& zgASX2)}a=|*LB@7yRTe#XLeafOZhk!EU8!UyYO+<`p zb{XtEO6+8gE9cbm%ZjCYVz^S_omk7*OJOhBZGse4dw7n>9{CIIc~y4XcsM3~ZXz7J z@PYjhcSUxWv3$Buw!ChQufCcX#5oIfGuyO>&~u#e|3ew;jim zQocqMJGiAfNI_4!plej+d4kj71`fY&TZY-6*omwaMJDWrJS+NOob}a1-?^4$zA<#{ zvy)+WWlmnziJ2LT@xtVI@N!#reSS}&vhRg#n#MuL5}M!hw|c{QU1H`58B9(?#oh&j z#TNB^e&#<54%TE+%m*o_bY#>Cf5eKvdaWl?QdL$5nK&XD$}8J(;|mSUs~lrwbixjG z6iw+e3JzvUkt%+f9t(P|vdZWiN&GcVt1wWvl-b<`ZLq;MtkjLIxSl9ZEs+qLH3pfR zu__`br(0lsNg9+2%N>p8y)qB&h~+6(Hd(v0IXLH07yUgyhM)y`Y^$l2_bubiQn`^B2Qq5UhH#EC7d zWI3$?_pRc$O(z>OTW+<++Jq7`arC``VkwJ;g(?Op@}PJXsmfrEClyH+lD2b8gJ+F$ zE1n|L%?B5WYva;Xsic`+C2r7Q9@XZx`(^XNkZ$uu$wh7PtYe^5cLK4~TPkup_c~} zvwM-5dAV?yF;%CNjnw$tFoB;~KiB7~YeO1)^j!9L=s!)(%d3~#z+Lxxgec*;fo>IY zNBl7tpS4O;4yR`6%0!Zd=Jkz4tUPTS*`||{wUWT&EoAm0E;_n171HFhv;7QEJ7u=D z5#b{StYP`lws1X*luIp7PxrGC0X;TafqlQa)n;y?&4V$q$(s3^pINr6${AI2@QPu3 z(bJm*Y7w1TrtZu5)EWn3-lv5)$;Pgc@?>>o-1fw>(FWpFFF_LPZqKpk)}#_6YW^Xyx~lBDxw&p|c8u8( zvNrc#d9I#K^yc2MRQnTLWZY1AVy-8KEY)pevpVHqD{x~p1{R;V&_h>Bm0n8T*C9)S zQ5zrUlog*iz5p*=S`@cPwK1rP9AB^#vwu5;9Lpo0b3f~h-E1^`$skAytF9Czswz6f zv$A3sD(`2dN}01^fH*^En@sd0aoP-ZZm6b z!qF37Et~KVvGZTyiCbLJSEuulHgSqQ6dwM}LD590*kP#9Mi%O4yZ$xK zMX6CFNylM<<#6*zbhvkfYJIOOL9$P%htlzlu*dCmX!UNhY;CH}A}pnDJ!&menWcu< z^?>oLdR;51N5odCk+C{PCzh(uc7alI%6VusUcK3F(WbhV%tBU256YY=Qk&}L2)ui9 ztBS43$`jlWso38xDhqD1>Q{!P6o_7UGs={$7k2s}iI_hMt&-pz(AU3P0n0YzDIMUOeYB~&aHMNFvXJE~qVQ-PCqe3tu z(86k%bCqQTsxu~CLy$5zlP`DlLic3;?7j1Bg}{*ScF5vyY^JxN)dcclwL!;mHO;nQ zc|*5`CazGcF~b^9PcAQE$F=PA((FquMA6|7@cc4#GC9sdM{nxDt+7aPM zD!;N}Oo)<>r~hQpIn2=3tujP#7|h_YXGp#Dc6!|Hq&Un0J^y6dZIix~%<{vj=-UxP z7Z?|>UQM^o=+aXzaV{>w_4OOEi#1;`D)lXNt9SUi$+ci}_5488G*Y0ac!fESg?6*{ zZIpsmY0`B0d3B-3&8moPENQdde&f*MU+ShEs%vH72BE>8#e z_ts`8qt3@(?$Qr5I-9R$EuOM)72A=ey?sp6NP|if6oZ9^3==0R#3vwek&vPnZ0KpJ zlH71lmZn-~v}T8H9@*==%Qy|W3d;+o8p+AFG5V+ldjEVvBLS*x4#xb7SK({#_;%r* zsDxq|v$X@pO^;>LLT$8&6iL00^`bU85KeGA*qaPbH{I)qPZ@3wwGwfaQV9RHRPrsy z1};z`&fV>K>!4JEIEpPLx+rulbie9Pm5uDFbzKUUsmnh&;tO}SeluMcS#G6+dm*wn z;vgBc>^fLXnR+ruZPB={eEUqa?edI}J-kW5)AA}ZEn$H1$T39U-OWUpciAKUdCHL_ z<(pN8p`HZJF$o(=ku>e2BKsbaku-gr_v`$RUKE9MYei$D@lJlLR4C)ysz8u22vKNm zy3@yy)$m!RxJ?!mA_Sg$5Cyn&GLT3Mrn|zyAOuZa_~*1y%=SlKlN^Uf5*mV%rz~YW2L1_R6z#2tJxA? zOCs^yml~*Y8cYfL8w1e;O+*fxbCZ^H@88-+Q7f+zk;JTzzp;-d;Xil#1U(LNBBDJO z^j>zn)r&YN^0N(r(v=%57f|icPV6#icLldOz*8BVMEi#EX>MYQk3CW^>#Oipog{g% z`mX5>o;S_A&Kpv6ud9TOpN00W)0+;PvGq!T^x2BuCm$!Tk<0z)8eEb8$onAxuVj3C zK4h!HDTM+a5C~tZ>|w!7Iz8KprrVonD4U2kc(a_Z1;+pjs9;(O&7Q4ox=zW4%BYMI zIqReGM9`-Y6N~Sze)TiPW0X$TaIGa^ih6{flgsOj#XR5gL>Bkt(;=Q}cEZ;b1TU;?2KJe%-Md=marmwHI^m+zgVO@6$ zGp#UMX`kFn!HjvI%zk{rdsv5k_%?cuE#I9w)!><6>iO>}-2GRW(i`{;o-A#GgBn5c zG9wC_T|WMw5;-`3T_l>wyiBBi_vom@lgFUqy7tHD)<>EKm!+V_xPz@yw{C9M-zUFJ z<$N)}=T%E42DH$Vq@#!9&gkTDLca^w>hHPboXZRbr@~X-Z#G;W5I3cAh#7YN%3qn) zYV%~7YH1JWcsXcM{2x0zYTk|M!)aIhI;8u_9*ocRD>~mt8-K356$p*?7 zy40|l^YuY43m@m+koSWA0psa zPnKF)700StS8=erByKKjAT)K}LZLRbOx;o0^LN=4ueuecT!fpD=G8E_-3JkaozBpz z;9{#zdqdpVy$0)_<>e+EY*KB9p33lGef?ea&F_wvLBABAL)Ulv?70X~DypH+CgX-= z!>t9QCj~y(>pDoaW=*1k!Y2l)B za^$)2s^M0UuQf|s<4`ugM1?^+PmJQa59nJ%P|9584frI4&++hPum{EWqY&dqE;L9` zGCH@>`+|#Cyo9jMF=2$9PS9;M~C`dEkTZ^VH)b4LR>!y1AQ=AFPF(ZHr($K+{7!OR-Pkx-j zcc6HFlb6!d7#L_QxpB_DRCCm}nMm?8)Wj+F&sDVbuX*Ev*gp74RoH>LHHFQ{>p(yw=Gh7KaQFGxd3gN7YjNEA0Rwp>r1a`BMnRpeC` z_awMGtTojUkcFuByaVa_JIA@L&FhAKWKqX-zr5qzjQ#b?VeKoEaBB|Q0`^vLsk>MT z6#BvUGSP6n&B1E#1fV9!^zAzIU&&?-Z}FxKAJ=OdB9W~=8XPrC`sM(!VDAD~>qSgBDw7GMq!q06rCL`bkMWflxod z1Rw_9^BX-vW&m&Ffmib$Z2;5h@d6ep8yh*#dpc5+3T}(B{5mkS%~h@iO1#2)6yy3P zDI3=q#iiv`v`%c`aWxLM3{;r<=Ppx3y+Y6u>ns`=balm;bvN8?>HhpR{Q@UocYcbpa zloXDlAc8$>dr7ZRYJ0u!P?w+&KX^aEyET%AODXtirS%&@$RS(=pq)6?7p5(3@d6^| zqOK%THLfV|`t|YP$OX!c1efb~fK@g_oL0^j^ywOL_+R#V9;Im%u`~8Xm)&cg41Kx|hvC^x_X`@2u z?=Pu=hc4dOE#>^h)7E+P)9bqe2K%^5PiE}wV7zN-p+Sxh0qW3s{eL{V009%TE7+Yk~$K*rx2J3VN1F+j)>#p#RIV^oX?_+3&rvgeD{ zf{BmOBHbe8iFkU7{^jxV)neMJFLj;!Ys?&KeA@jo>%)ht-%1l(tIAFTogcce&?1>u z`ILslb?=MQZI!Cba5Efe&T0bw$a44l))(dH;@z#&`05@W9)`bs)4B^j8*RuR!?e0Cq`+SBtfqtE?w^6@kM} z$wx%5Hs`F7rJX2TGrfjB1o*&}(I1XHmag#vl;@yvI;H?5KC%L1_UpH$9&o93rpJX5 z9|9kI=J``<2=q8n8Vm`fXj+uz60n^oB@LzE^JL+1!=X+ObaeBP)`RW>m&| zyd+I3%$ymlx_q|P8h}fDMH5i}LKHPy{27J8eR9PBd!2w)ZrTSWz8@B(_PbbS@xenU zLf_?{zQ*w$14+C+mU#AQd-2LB!(?AX z_GSk0=Vu}S?qmCN566S6F{26l$}n+!2{hXCwuEPK(m@s4JZ}8r>C`vxK$7z)Odwg8 z;C9GjVXdTA^jx}QzQ@I3@~E0t_}cGxRu2?Pv7;C9zp0PlNo&{l;$xA*Nx%D%<&Fh; z2|!E5J#TK7zdGC!^t||v>2Y8TYB2{ev%kLwgP)g5?{*QZLNc1HWBB_BUt!YsPTw%| z%z-G^3)FHHCi9g0hEs$c{SSWpEI(+#&ZoH+)=#OISzL^uU9mx6#wOz zUjlsn*Q_3pI*SJm>yKM8^XFGb0I3%F%Gg;)&QZ>b%0iGfF~E#}_Q0jFF28UM0+Ykv zjo{Q`X;|{ON=?h#x&5^E=N(b&(eHbv0O|>I3lhx-!j}eb*Kwr4B$1N}!nN22Hdp#U zDSY`Yj!rLo{K^s$_~0t^E}Ur}h+BZK-~mWq4JSWlZ(ivG_IZWB*uIVi6-Da(2^waX zym7nncXm=nBgevv!y~F;Eh-XaEN11>0h@*cgB-?^4UDP1c(}-&>4ttTJCg*prhXt- zi24-Pp1H_LCWRRfLG49QJ&oA@*P!M?c~;N=kzB!+w3z5^WNiQ?dh<={|w z^ZMVJBTtT+@dq!_Ho}a8iMkHF$M=35g(-IR`<*8(Y}%s&4S>E0K*p5eX&Lf-7peam zyjKETE*35Dh3ve})fM))?@39+x{@hiEI$gmdcDH~WODTG@8D0Mp-ALNg25MohCTmK z>P#@$1t6&bj73$g2Rieo#tXaq6*MqPn*2?{YtvN6()VPY1iJ;7I#KQ!yQHq*fS9Wt zj}wopmso!k6#5dhzst~(#CmHi^<&$w=40BMbOn>LYiD{aEoUv^b&p<vVJ3>Ye4u-xvQ|LAN5R5k2H5)crX5rm|=ZKavG0?lDWeny~p z({8vyFx@0jH&7Vj07kDtB6DUCb8s`z+wA#)gkZX7TQXZhv-M~HIaYY{@6wZ{G~A?X zOY6@(F_IKCIr!@GAWKAyqykMSz8V5$wNp_$ZVvdl9o5sPIHn55IowLz96U-7TW+P7 zz_*uS?4*Ey8espCk&zLC@lx$Ooub5^B!r!RHSq+NHA1q@aNgOB9gxT+5jL3WZvvYr znD%f_aq|+rpj>})2jS67FzK%-B#G*bDaDCxS6D{Zhm>GWWAK-$>k$GqFUsPbx1Y9X zu`2zE{zAOeS=GhYR@KM%u8xH&Jmzp?YXw&zKTWPR3N#KER{Q};1c;(3oH17X-$wO;725}!4jMUUGuDt%bwubRbJjxH$ zAV6YkVzB!1Lv&A+T;rP5ofsVcM*xQCq9%oDqKU3aBJijCVJaqGXNpYQ!1UfY5B&UB z?j=KWkiD_t0Bvp8gY+IbH8w9)W&}EOEBmwj?5r1J5+UAP)!$c7@o}n98hq=-2xtPd zGQ4pG?xwQP#q7%mQ2@q72EhVBaV{rB&`B>tz)hG~fD64}#IG16YkS)%p9^@~$p3(H z`ND93{Q)=-i;BMEQ`5dEYW|wJ%X9*KFwFxGA0`2UuzJcf*^yz3D;BO0Eo0x7yvfT6 zo6~J9wrk{SK{go}Sim(e4sado))3ptjmmpbySoNcjSz7CaPge8o!!RU26@N5Acgod zyWQc2WCOg?ExKI6czRBK)%~C#in0M>0bhR|1K7H$;kaS5`^CO6AR-)bD#*eHHRgO7 zYy$kq_)-)JLmFvU!_70dav+#qa;l}=5+*Z%jXq4cLh7rtf6n@!9_C4lSIx&212}_& z*WKpb!e_}_!Y?eKxj>lxg&2V%?Cpbwo-h0rBK%bXhKXf5n;9C!oHhAv+5muLDSNIt zKpkBb`U%*$|K_|H5C~jbj%N8h*cZ(ZAS(F7bUzIEH`a%$fuQl>>^k^?*e^}|SyB=m zf$bo<=>7cs9liHa7ngW^ckuq|_1KlEtcii|Hm(`{>+tWU+&KQK!HnFIta&oJ|0uW= ze3o6oWX})?G}>TN&H^-y=BBBoR-nRVFS4G)1=x_L1~9kPZ-~d=34ixsEBUACM`|FX zJo)gyXL{miR@bWWr@u=YFG(~twX8H#Brp*HSsM`;8@?DDK9!n=i^#BigIl;)9_hnO zKSd#>6nqcxe-m{8d$s}C_-#bD=6gV|Vxu>6ih}@E-5_TCojAa(?W##g0Bq>cgZj;P z1MKAwhKC>w2qlC`bxfR4(LfL>z|oj#c_Up<3_(FnyG$suu7Ne0)Jk*Mw?c;J8>8w! zNWPmI^S&e}ee~=r6=z(!=GcyNg8^sp+r07ZQ^2emUUiTeH~}gX5X1soO)WZjBlMw< zkwyQu8TN0MzCU~`{E9#&dALIJ3G}iCz^Lw2wRS(?2=+?f53uuYMdH)8K%RIa{IAgj zAA;;PTlIeq?H7=#ITGacj0O!%uPXEq^yH-8I|JD%Q@rL6HnR)JK%g9v1N8pY2L$E0 zTO9-Jn~?+<4>i5_#Cj2R{WLsWtYFf&;hmH(591?ieOXq{cdUTtlxrdB%U zqv~>1nGIqmO8|>bfKLC##~%MTKAr$fcdDx%G}o)i2Opx$1E^uz91Ym}*TR6_yt`z7hJkyc$`+~l?lzBO`m^(Y&$q&=qNRYMlC! z*L`kd;)9I)s!M3+DnSW|q!pYhEslZrK8};EL%!m(u{!3aW;oHTrez=r1dev?rxmAq z4`rteXiF};c4-8Su{TKKfv{bB1gyFo^D{~3Y*Yy-P*8ydB;ba!eq^odm<1v8f!r!g z{3paZJkrsO??_2?&IrH|2(LIU?o;J>5TFZQ6e>PQA`)vWl$7??LfV=$e2N?1F9Ej&Urokj2EM@i~0*0v=frjEkU?!>Ffq=h5JL923 zV|fOt$2X<6dxAx!&e@8+w^x^^l`a_FF!{b@f4eTM(WxS;M5NKs&5fvm6}Hd?bAegC zE2-6D86B&ijxI>a&x4aqJF2SRJgFCgTW1ZVU!}860r??1X5XFp_}Kb(h#|@}tX*o{ zh^}Yo_$FV`o^!+vm={#e087PEQ6K^0_z#epCiAb?BnZ3y878TWrOm+yd#;P@=JWYy ztni1{{wq>Mv(Y0;-ujX<>~FhhZVqJ0z+?bByfMuqX5 z`oOaLhd&QGE~N;;f47<)h$X=BLG6M7k1t{W=t9rqcDPjqZXR12E;Zo1X&nqh_Vm z>+ZI$!xsZv?2qTRbALRLK5WgU@HQC~~IKe54%k_D5b< z{l7zjF`5h1IBz8*xFGZz%pD64d{Qfa2|38Bp{j5piL)kfY-gi#4b0Od__9r)ohcDX zo`9p{i~i-ymo}2t*4E}))#g)i+1``a&C)OYVAJ;Yw!SE_nlKcQ0htHsCAR(!#rGf!W{iGj>!1Jk+UTi>#a_Li7ejwvX)G*BRW{7})lz;$yaQWL) zdC89aY-`xu{4@7B#yE{su0(+1*C@&;+x@kNY?)^(Y2u!zyAN3=AnfMe3i->sXgF3t;C@z%z`Dxf%^ zdNh!0G$P7#edGA=xbnLc%acw3WG+p-=4ZV8#Wp7Pfd}4~ZQj(s5O}syqrk$H4p_JX z1HQX4@DOaP15c{{O`sR{m5WE@Y50SKL(xinbCpWC=LTy){2}n18|bbVWnXzV))IC; zx1kach;=!wVdG@A0@H@g4MfRo22@m!r&RNiTQ5S5H zs3T*#;9zqi=d}{+_pg%PWioW&Zsv6Vo{xpyjV=dswVyVD5ORuLyQ&=x1s_;oFjkjK zzw&vobiC$64;3JheOt#s{I~){n*9%lMIV7E;K<*01~9k;1bqJvJZ}KMdEITX-}Zc& z;OVi1{k|-TQ1>J*VVlIvTc#r>Q-}BkZf3kMu_WAIN= candidate_start & solar_time <= modeled_end, ] + max_lag <- max(round(v(candidate$travel.time) / timestep_days)) + lead_in_start <- modeled_start - as.difftime(max_lag * timestep_days, units='days') + + full_data[solar_time >= lead_in_start & solar_time <= modeled_end, ] +} + +test_that("metab() fits a two-station model and predict_metab()/predict_DO() work", { + skip_on_cran() + skip_if_not_installed('rstan') + + small_dat <- subset_2station_data(two_station_example, n_modeled_days=3) + + sp <- specs( + mm_name('bayes_2station'), + n_chains=1, n_cores=1, burnin_steps=100, saved_steps=100, verbose=FALSE) + + mm <- metab(specs=sp, data=small_dat) + expect_s4_class(mm, 'metab_2station') + + pm <- predict_metab(mm) + expect_s3_class(pm, 'data.frame') + expect_true(all(c('GPP','ER','K600') %in% names(pm))) + expect_equal(nrow(pm), 3) + + pdo <- predict_DO(mm) + expect_s3_class(pdo, 'data.frame') + expect_true(all(c('DO.obs.down','DO.mod.down') %in% names(pdo))) }) From c76e7d6130a9814eec37d67c3b31f7b8af5e723a Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 13 Jul 2026 12:11:13 -0700 Subject: [PATCH 07/17] Fix R CMD CHECK NOTEs in metab_2station.R - Add utils::globalVariables for ., metab_50pct, DO.mod.down - Add @importFrom utils modifyList - Replace require(rstan) with requireNamespace() pattern --- R/metab_2station.R | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/R/metab_2station.R b/R/metab_2station.R index e298dd5..6727472 100644 --- a/R/metab_2station.R +++ b/R/metab_2station.R @@ -1,6 +1,11 @@ #' @include metab_model-class.R metab_bayes.R NULL +# Suppress R CMD CHECK NOTEs for column names used as unbound globals in +# dplyr NSE calls (e.g. mutate, summarise). These are data frame column +# names resolved at runtime, not missing variable declarations. +utils::globalVariables(c(".", "metab_50pct", "DO.mod.down")) + #' Two-station Bayesian metabolism model fitting function #' #' Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to @@ -33,6 +38,7 @@ NULL #' #' @export #' @family metab_model +#' @importFrom utils modifyList metab_2station <- function( specs=specs(mm_name('bayes_2station')), data=mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, @@ -126,9 +132,7 @@ metab_2station <- function( inst <- NULL withCallingHandlers( tryCatch({ - if(!suppressPackageStartupMessages(require(rstan))) { - stop("the rstan package is required for Stan MCMC models") - } + if (!requireNamespace("rstan", quietly = TRUE)) stop("rstan is required but not installed. Install it with: install.packages('rstan')") consolelog <- utils::capture.output( stanfit <- rstan::stan( From a6eb23b20f3590d62342bc52dbc55ac28b4aa43d Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 27 Jul 2026 12:49:13 -0700 Subject: [PATCH 08/17] Rename metab_2station -> metab_bayes_2s, mm_ts_prep_data -> prepdata_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_()-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). --- DESCRIPTION | 5 +- NAMESPACE | 12 +- R/data.R | 18 +- R/metab.R | 2 +- R/{metab_2station.R => metab_bayes_2s.R} | 249 ++++++++++++------ R/mm_name.R | 60 ++--- R/mm_parse_name.R | 47 ++-- R/mm_ts_prep_data.R | 125 --------- R/mm_valid_names.R | 4 +- data-raw/two_station_example.R | 17 +- man/get_params.Rd | 41 ++- man/metab_Kmodel-class.Rd | 1 + man/metab_Kmodel.Rd | 1 + man/metab_bayes-class.Rd | 1 + man/metab_bayes.Rd | 1 + ...ation-class.Rd => metab_bayes_2s-class.Rd} | 11 +- man/{metab_2station.Rd => metab_bayes_2s.Rd} | 37 ++- man/metab_mle-class.Rd | 1 + man/metab_mle.Rd | 1 + man/metab_model-class.Rd | 1 + man/metab_night-class.Rd | 1 + man/metab_night.Rd | 1 + man/metab_sim-class.Rd | 1 + man/metab_sim.Rd | 1 + man/mm_generate_mcmc_file.Rd | 17 +- man/mm_name.Rd | 19 +- man/mm_valid_names.Rd | 19 +- man/predict_DO.Rd | 30 ++- man/predict_metab.Rd | 18 +- ...m_ts_prep_data.Rd => prepdata_bayes_2s.Rd} | 28 +- man/two_station_example.Rd | 18 +- ...metab_2station.R => test-metab_bayes_2s.R} | 56 ++-- 32 files changed, 447 insertions(+), 397 deletions(-) rename R/{metab_2station.R => metab_bayes_2s.R} (60%) delete mode 100644 R/mm_ts_prep_data.R rename man/{metab_2station-class.Rd => metab_bayes_2s-class.Rd} (75%) rename man/{metab_2station.Rd => metab_bayes_2s.Rd} (65%) rename man/{mm_ts_prep_data.Rd => prepdata_bayes_2s.Rd} (67%) rename tests/testthat/{test-metab_2station.R => test-metab_bayes_2s.R} (82%) diff --git a/DESCRIPTION b/DESCRIPTION index d3852d2..5607700 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -95,9 +95,9 @@ Collate: 'metab_model_interface.R' 'specs-class.R' 'metab_model-class.R' - 'metab_bayes.R' - 'metab_2station.R' 'metab_Kmodel.R' + 'metab_bayes.R' + 'metab_bayes_2s.R' 'metab_inputs.R' 'metab_mle.R' 'metab_model.get_param_names.R' @@ -125,7 +125,6 @@ Collate: 'mm_predict_DO_1ply.R' 'mm_predict_metab_1ply.R' 'mm_sd_to_ci.R' - 'mm_ts_prep_data.R' 'mm_validate_data.R' 'plot_DO_preds.R' 'plot_distribs.R' diff --git a/NAMESPACE b/NAMESPACE index fc43e76..eb89725 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,21 +10,21 @@ S3method(get_mcmc,metab_bayes) S3method(get_mcmc_data,metab_bayes) S3method(get_param_names,character) S3method(get_param_names,metab_model) -S3method(get_params,metab_2station) S3method(get_params,metab_Kmodel) S3method(get_params,metab_bayes) +S3method(get_params,metab_bayes_2s) S3method(get_params,metab_model) S3method(get_params,metab_sim) S3method(get_specs,metab_model) S3method(get_version,metab_model) -S3method(predict_DO,metab_2station) S3method(predict_DO,metab_Kmodel) +S3method(predict_DO,metab_bayes_2s) S3method(predict_DO,metab_model) S3method(predict_DO,metab_night) S3method(predict_DO,metab_sim) -S3method(predict_metab,metab_2station) S3method(predict_metab,metab_Kmodel) S3method(predict_metab,metab_bayes) +S3method(predict_metab,metab_bayes_2s) S3method(predict_metab,metab_model) S3method(print,logs_metab) S3method(print,specs) @@ -70,9 +70,9 @@ export(lookup_google_timezone) export(lookup_timezone) export(lookup_usgs_elevation) export(metab) -export(metab_2station) export(metab_Kmodel) export(metab_bayes) +export(metab_bayes_2s) export(metab_inputs) export(metab_mle) export(metab_model) @@ -86,7 +86,6 @@ export(mm_model_by_ply) export(mm_model_by_ply_prototype) export(mm_name) export(mm_parse_name) -export(mm_ts_prep_data) export(mm_valid_names) export(mm_validate_data) export(mm_validate_name) @@ -99,9 +98,9 @@ export(revise) export(sim_Kb) export(sim_pred_Kb) export(specs) -exportClasses(metab_2station) exportClasses(metab_Kmodel) exportClasses(metab_bayes) +exportClasses(metab_bayes_2s) exportClasses(metab_mle) exportClasses(metab_model) exportClasses(metab_night) @@ -170,6 +169,7 @@ importFrom(utils,available.packages) importFrom(utils,capture.output) importFrom(utils,contrib.url) importFrom(utils,head) +importFrom(utils,modifyList) importFrom(utils,packageVersion) importFrom(utils,read.csv) importFrom(utils,tail) diff --git a/R/data.R b/R/data.R index cc45a05..b0a00ab 100644 --- a/R/data.R +++ b/R/data.R @@ -1,17 +1,17 @@ #' Example two-station (VFTS) input data #' #' A 30-day example dataset for fitting a two-station (upstream/downstream, -#' a.k.a. VFTS) metabolism model with \code{\link{metab_2station}}. It is a -#' subset of the \code{VFTS-2} (variable-travel-time) run from the published -#' two-station metabolism modeling dataset for a reach of the Colorado River -#' in Glen Canyon, covering 2011-07-31 through 2011-08-29 at the source -#' data's native 15-minute timestep, plus a lead-in block of upstream DO -#' observations (exactly as long as the longest travel time in the dataset -#' requires -- see \code{\link{metab_2station}}'s "Two-station data -#' requirements" section) immediately before 2011-07-31. +#' Variable Flow Two-Station) metabolism model with +#' \code{\link{metab_bayes_2s}}. It is a subset of the \code{VFTS-2} +#' (variable-travel-time) run from the published two-station metabolism modeling +#' dataset for a reach of the Colorado River in Glen Canyon, covering 2011-07-31 +#' through 2011-08-29 at the source data's native 15-minute timestep, plus a +#' lead-in block of upstream DO observations (exactly as long as the longest +#' travel time in the dataset requires -- see \code{\link{metab_bayes_2s}}'s +#' "Two-station data requirements" section) immediately before 2011-07-31. #' #' @format A data.frame with 2904 rows and the 9 columns expected by -#' \code{\link{metab_2station}}'s \code{data} argument, each carrying +#' \code{\link{metab_bayes_2s}}'s \code{data} argument, each carrying #' \code{\link[unitted]{unitted}} units matching \code{\link{mm_data}}: #' \describe{ #' \item{solar.time}{POSIXct timestamp, UTC} diff --git a/R/metab.R b/R/metab.R index 6e6b705..4d55efb 100644 --- a/R/metab.R +++ b/R/metab.R @@ -60,7 +60,7 @@ metab <- function(specs=specs(mm_name()), data=v(mm_data(NULL)), data_daily=v(mm metab_fun <- switch( model_type, bayes = metab_bayes, - bayes_2station = metab_2station, + bayes_2s = metab_bayes_2s, Kmodel = metab_Kmodel, mle = metab_mle, night = metab_night, diff --git a/R/metab_2station.R b/R/metab_bayes_2s.R similarity index 60% rename from R/metab_2station.R rename to R/metab_bayes_2s.R index 6727472..8b192d8 100644 --- a/R/metab_2station.R +++ b/R/metab_bayes_2s.R @@ -8,22 +8,22 @@ utils::globalVariables(c(".", "metab_50pct", "DO.mod.down")) #' Two-station Bayesian metabolism model fitting function #' -#' Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to -#' estimate GPP, ER, and K600 from paired upstream and downstream DO, -#' temperature, light, and travel-time data, using the single fixed Stan -#' model in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See -#' \code{\link{mm_name}} to choose a Bayesian model and \code{\link{specs}} -#' for relevant options for the \code{specs} argument. +#' Fits a two-station (upstream/downstream, Variable Flow Two-Station) Bayesian +#' model to estimate GPP, ER, and K600 from paired upstream and downstream DO, +#' temperature, light, and travel-time data, using the single fixed Stan model +#' in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See \code{\link{mm_name}} to +#' choose a Bayesian model and \code{\link{specs}} for relevant options for the +#' \code{specs} argument. #' #' Unlike \code{\link{metab_bayes}}, which supports many model structures via -#' \code{split_dates}/\code{pool_K600}/etc., \code{metab_2station} always +#' \code{split_dates}/\code{pool_K600}/etc., \code{metab_bayes_2s} always #' fits every date jointly in a single Stan call (\code{specs$split_dates} is #' forced to \code{FALSE} by \code{\link{specs}}), because the #' upstream-downstream lag shift ties each date's first modeled rows to the #' previous date's last rows. #' #' @inheritParams metab -#' @return A metab_2station object containing the fitted model. This object +#' @return A metab_bayes_2s object containing the fitted model. This object #' can be inspected with the functions in the #' \code{\link{metab_model_interface}} and also \code{\link{get_mcmc}}. #' @@ -39,8 +39,8 @@ utils::globalVariables(c(".", "metab_50pct", "DO.mod.down")) #' @export #' @family metab_model #' @importFrom utils modifyList -metab_2station <- function( - specs=specs(mm_name('bayes_2station')), +metab_bayes_2s <- function( + specs=specs(mm_name('bayes_2s')), data=mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, depth, temp.water, travel.time), data_daily=mm_data(date, optional='all'), @@ -49,46 +49,24 @@ metab_2station <- function( stanfit <- NULL fitting_time <- system.time({ - # Check data for correct column names & units - dat_list <- mm_validate_data(data, data_daily, 'metab_2station') - data_v <- v(dat_list$data) + # Check data for correct column names, units, and travel.time bounds + # (mm_validate_data()), then check lead-in coverage + # (mm_validate_data_2station()), before any data prep begins + dat_list <- mm_validate_data(data, data_daily, 'metab_bayes_2s') + mm_validate_data_2station(dat_list$data) + data_v <- v(dat_list$data) travel_time <- data_v$travel.time solar_time <- data_v$solar.time - # a. travel.time must be strictly positive - if(any(travel_time <= 0)) { - stop('travel.time must be > 0') - } - - # b. travel.time is expected in days, so any value >= 1 almost certainly - # reflects a units mistake (e.g., minutes or hours rather than days) - if(any(travel_time >= 1)) { - stop('travel.time must be < 1 day; values >= 1 suggest incorrect units (expected days)') - } - - # c. there must be enough lead-in rows of upstream DO before the first - # modeled row to cover the longest travel time in the dataset. timestep_days - # is the median observation interval, in days; max_lag is the number of - # timesteps by which upstream data must lead downstream predictions. The - # first max_lag rows of data serve only as lead-in and cannot themselves be - # modeled, so at least max_lag + 1 rows are required overall. + # Reconstruct the same "modeled rows" (post-lead-in-trim) index set that + # prepdata_bayes_2s() computes internally, using the identical + # timestep_days/lag/max_lag formula, so that Stan's date_index/time_index + # can be mapped back to actual dates and solar.times for the daily and + # instantaneous results below. (Duplicated here rather than exposed by + # prepdata_bayes_2s(), which returns only the Stan-ready matrices.) timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) max_lag <- max(round(travel_time / timestep_days)) - if(nrow(dat_list$data) <= max_lag) { - lead_in_needed <- max_lag - nrow(dat_list$data) + 1 - stop(paste0( - 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', - max_lag, ' timestep(s), but only ', nrow(dat_list$data), ' row(s) were supplied; ', - 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) - } - - # Reconstruct the same "modeled rows" (post-lead-in-trim) index set that - # mm_ts_prep_data() computes internally, using the identical timestep_days/ - # lag/max_lag formula, so that Stan's date_index/time_index can be mapped - # back to actual dates and solar.times for the daily and instantaneous - # results below. (Duplicated here rather than exposed by - # mm_ts_prep_data(), which returns only the Stan-ready matrices.) n_total <- nrow(dat_list$data) keep <- seq.int(max_lag + 1, n_total) modeled_solar_time <- solar_time[keep] @@ -109,11 +87,11 @@ metab_2station <- function( time_index=rep(seq_len(n_obs), times=n_days)) # Prepare the Stan data list (matrices from data, plus scalar priors from - # specs). modifyList (not c()) is used because mm_ts_prep_data() already - # supplies K600_lnorm_meanlog/K600_lnorm_sdlog (with its own fallback - # defaults), and those two names are also in specs$params_in; a plain - # c() would create duplicate-named list elements instead of overriding - data_list <- mm_ts_prep_data(dat_list$data, specs=specs) + # specs). modifyList (not c()) is used because prepdata_bayes_2s() already + # supplies K600_lnorm_meanlog/K600_lnorm_sdlog (read from specs), and + # those two names are also in specs$params_in; a plain c() would create + # duplicate-named list elements instead of overriding + data_list <- prepdata_bayes_2s(dat_list$data, specs=specs) data_list <- modifyList(data_list, specs[specs$params_in]) # Check and parse model file path @@ -201,7 +179,7 @@ metab_2station <- function( # Package and return results mm <- metab_model( - model_class="metab_2station", + model_class="metab_bayes_2s", info=info, fit=fit, log=NULL, @@ -230,11 +208,133 @@ metab_2station <- function( } -#### metab_2station class #### +#' Reshape long-format two-station data into the list expected by the +#' two-station Stan model +#' +#' Time-shifts the upstream DO series to match the travel time between +#' stations, then pivots the result into the \code{n_obs x n_days} matrices +#' expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} +#' (see \code{\link{metab_bayes_2s}}). +#' +#' The upstream observation that "matches" a given downstream observation at +#' row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / +#' timestep_days)} timesteps earlier, where \code{timestep_days} is the +#' median timestep of \code{data$solar.time}. This must be computed the same +#' way as in \code{\link{mm_validate_data_2station}}'s lead-in check, so that the +#' \code{max(lag)} computed here always agrees with the lead-in requirement +#' already validated there. The first \code{max(lag)} rows of \code{data} are +#' lead-in rows: they supply upstream DO for the shift but are never +#' themselves treated as modeled (downstream) observations. +#' +#' @param data data.frame as validated by \code{\link{mm_validate_data}} for +#' \code{\link{metab_bayes_2s}}: must contain \code{solar.time}, +#' \code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, +#' \code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, +#' \code{travel.time}, sorted ascending by \code{solar.time}, and must +#' include the lead-in rows required to cover the longest travel time (see +#' \code{\link{metab_bayes_2s}}). +#' @param specs a list of model specs (see \code{\link{specs}}), expected to +#' already contain \code{K600_lnorm_meanlog} and \code{K600_lnorm_sdlog} +#' -- e.g., the object returned by \code{specs(mm_name('bayes_2s'))}, which +#' populates them with sensible defaults. This function does not supply +#' its own fallback values; if \code{specs} is omitted or missing these +#' fields, the resulting Stan data will contain NULL/missing values for +#' them. +#' @return a named list with all variables in the Stan model's data block: +#' \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, +#' \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, +#' \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} +#' matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} +#' @importFrom unitted v +#' @keywords internal +prepdata_bayes_2s <- function(data, specs=NULL) { + + # strip units; Stan cannot handle unitted vectors/matrices + data <- v(data) + + # timestep_days must match mm_validate_data_2station()'s lead-in check + # exactly (median timestep, in days), so that max_lag here agrees with + # what was validated there + timestep_days <- stats::median(as.numeric(diff(data$solar.time), units='days')) + + # lag, in timesteps, that the upstream series must be shifted by to line up + # with each row's downstream observation + lag <- round(data$travel.time / timestep_days) + max_lag <- max(lag) + n_total <- nrow(data) + + # the first max_lag rows are lead-in only (upstream DO used for the shift, + # but never modeled themselves). every row i > max_lag is guaranteed to + # have a valid shift target (i - lag[i] >= 1) because lag[i] <= max_lag + keep <- seq.int(max_lag + 1, n_total) + shift_idx <- keep - lag[keep] + if(any(shift_idx < 1)) { + # should be unreachable given max_lag's definition; guards against + # programming errors rather than expected user input + stop('internal error: upstream shift index falls before the first row of data') + } + + modeled <- data.frame( + solar.time = data$solar.time[keep], + DO_obs_up = data$DO.obs.up[shift_idx], + DO_sat_up = data$DO.sat.up[shift_idx], + DO_obs_down = data$DO.obs.down[keep], + DO_sat_down = data$DO.sat.down[keep], + light = data$light[keep], + depth = data$depth[keep], + temp_water = data$temp.water[keep], + travel_time = data$travel.time[keep] + ) + + # pivot into n_obs x n_days matrices, one column per unique date, following + # the same time_by_date_matrix approach used for the 1-station Stan models + # (see prepdata_bayes() in metab_bayes.R) + date_vec <- as.character(as.Date(modeled$solar.time)) + date_table <- table(date_vec) + n_days <- length(date_table) + n_obs_per_day <- unique(unname(date_table)) + if(length(n_obs_per_day) > 1) { + stop( + 'dates have differing numbers of modeled rows after lead-in removal; ', + 'observations cannot be combined into a matrix: ', + paste(sprintf('%s (%d rows)', names(date_table), date_table), collapse=', ')) + } + n_obs <- n_obs_per_day + + to_matrix <- function(vec) matrix(vec, nrow=n_obs, ncol=n_days, byrow=FALSE) + + # confirm each date occupies a contiguous block of rows, i.e., that data + # was sorted by solar.time; otherwise the matrix pivot below would silently + # scramble which rows belong to which date + date_mat <- to_matrix(date_vec) + unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) + if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { + stop('data must be sorted by solar.time so that each date occupies a contiguous block of rows') + } + + list( + n_obs = n_obs, + n_days = n_days, + DO_obs_up = to_matrix(modeled$DO_obs_up), + DO_sat_up = to_matrix(modeled$DO_sat_up), + DO_obs_down = to_matrix(modeled$DO_obs_down), + DO_sat_down = to_matrix(modeled$DO_sat_down), + light = to_matrix(modeled$light), + depth = to_matrix(modeled$depth), + temp_water = to_matrix(modeled$temp_water), + travel_time = to_matrix(modeled$travel_time), + K600_lnorm_meanlog = specs$K600_lnorm_meanlog, + K600_lnorm_sdlog = specs$K600_lnorm_sdlog + ) +} + + +#### metab_bayes_2s class #### -#' Metabolism model fitted by two-station (VFTS) Bayesian MCMC +#' Metabolism model fitted by two-station (Variable Flow Two-Station) Bayesian +#' MCMC #' -#' \code{metab_2station} models use Bayesian MCMC methods to fit values of +#' \code{metab_bayes_2s} models use Bayesian MCMC methods to fit values of #' GPP, ER, and K600 from paired upstream/downstream DO curves. This class #' inherits from \code{metab_bayes} (same \code{log}/\code{mcmc}/ #' \code{mcmc_data}/\code{compile_time} slots, and therefore the same @@ -243,26 +343,25 @@ metab_2station <- function( #' \code{predict_DO} are overridden below because two-station's fitted-value #' structure and output columns differ from one-station's. #' -#' @exportClass metab_2station +#' @exportClass metab_bayes_2s #' @family metab.model.classes -setClass("metab_2station", contains="metab_bayes") +setClass("metab_bayes_2s", contains="metab_bayes") #' @describeIn get_params Does the same Stan-output-to-streamMetabolizer #' renaming as \code{get_params.metab_bayes}, but (unlike that method) does #' not delegate the rest of the work to \code{get_params.metab_model} via #' \code{NextMethod()}: that generic implementation looks up parameter -#' names via \code{get_param_names()}, which assumes a -#' \code{metab_()}-named model constructor (would look for -#' \code{metab_bayes_2station}, which doesn't exist -- our constructor is -#' \code{metab_2station}) and streamMetabolizer's ODE-based dDOdt framework -#' for one-station models, neither of which apply to the two-station -#' steady-state model. \code{fixed} column/star annotations (relevant only -#' to models that can take fixed daily parameters from \code{data_daily}) -#' are not supported here. +#' names via \code{get_param_names()}, which builds an ODE-based dDOdt +#' function using streamMetabolizer's one-station instantaneous-rate +#' framework (\code{ode_method}/\code{GPP_fun}/\code{ER_fun}/ +#' \code{deficit_src}) -- machinery that doesn't apply to the two-station +#' steady-state model's daily GPP/ER/K600 parameters. \code{fixed} +#' column/star annotations (relevant only to models that can take fixed +#' daily parameters from \code{data_daily}) are not supported here. #' @export #' @import dplyr -get_params.metab_2station <- function( +get_params.metab_bayes_2s <- function( metab_model, date_start=NA, date_end=NA, uncertainty=c('sd','ci','none'), messages=TRUE, ...) { uncertainty <- match.arg(uncertainty) @@ -314,7 +413,7 @@ get_params.metab_2station <- function( #' the two-station Stan model results. #' @export #' @import dplyr -predict_metab.metab_2station <- function(metab_model, date_start=NA, date_end=NA, ...) { +predict_metab.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ...) { Var1 <- Var2 <- '.dplyr.var' fit.names <- expand.grid(c('50pct','2.5pct','97.5pct'), c('GPP_daily','ER_daily','K600_daily'), stringsAsFactors=FALSE) %>% @@ -352,20 +451,20 @@ predict_metab.metab_2station <- function(metab_model, date_start=NA, date_end=NA } -#' @describeIn predict_DO Two-station (VFTS) models. Returns a data.frame -#' with columns \code{solar.time}, \code{DO.obs.down} (the observed -#' downstream DO from the input data), and \code{DO.mod.down} (the -#' posterior median of the two-station Stan model's fitted downstream DO) +#' @describeIn predict_DO Two-station (Variable Flow Two-Station) models. +#' Returns a data.frame with columns \code{solar.time}, \code{DO.obs.down} +#' (the observed downstream DO from the input data), and \code{DO.mod.down} +#' (the posterior median of the two-station Stan model's fitted downstream DO) #' -- unlike the one-station \code{predict_DO} methods, which return -#' \code{DO.obs}/\code{DO.mod}. The values are those computed once at -#' fitting time (see \code{\link{metab_2station}}); \code{use_saved=FALSE} -#' (on-demand recomputation from the fitted daily GPP/ER/K600 medians) is -#' not implemented. +#' \code{DO.obs}/\code{DO.mod}. The values are those computed once at fitting +#' time (see \code{\link{metab_bayes_2s}}); \code{use_saved=FALSE} (on-demand +#' recomputation from the fitted daily GPP/ER/K600 medians) is not +#' implemented. #' @export -predict_DO.metab_2station <- function(metab_model, date_start=NA, date_end=NA, ..., use_saved=TRUE) { +predict_DO.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ..., use_saved=TRUE) { if(!isTRUE(use_saved)) { - stop("predict_DO(use_saved=FALSE) is not implemented for metab_2station; only the fitted-time DO.mod.down values are available") + stop("predict_DO(use_saved=FALSE) is not implemented for metab_bayes_2s; only the fitted-time DO.mod.down values are available") } inst <- metab_model@fit$inst diff --git a/R/mm_name.R b/R/mm_name.R index 4df02fc..17079dd 100644 --- a/R/mm_name.R +++ b/R/mm_name.R @@ -54,16 +54,16 @@ #' @param type character. The model type. Options: \itemize{ \item \code{mle}: #' maximum likelihood estimation (see also \code{\link{metab_mle}}) \item #' \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -#' \code{bayes_2station}: two-station (upstream/downstream, a.k.a. VFTS) -#' Bayesian model with a single fixed structure (see also -#' \code{\link{metab_2station}}) \item -#' \code{night}: nighttime regression (see also \code{\link{metab_night}}) -#' \item \code{Kmodel}: regression of \emph{daily} estimates of -#' \code{K600.daily} versus discharge, time, etc., usually for 3-phase -#' estimation of K alone (by MLE or nighttime regression), K vs discharge -#' (using this model), and then GPP and ER with fixed K (by MLE) (see also -#' \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -#' 'data' for testing other models (see also \code{\link{metab_sim}}) } +#' \code{bayes_2s}: two-station (upstream/downstream, Variable Flow +#' Two-Station) Bayesian model with a single fixed structure (see also +#' \code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +#' also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +#' \emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +#' usually for 3-phase estimation of K alone (by MLE or nighttime regression), +#' K vs discharge (using this model), and then GPP and ER with fixed K (by +#' MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +#' \code{DO.obs} 'data' for testing other models (see also +#' \code{\link{metab_sim}}) } #' @param pool_K600 character. [How] should the model pool information among #' days to get more consistent daily estimates for K600? Options (see Details #' for more): \itemize{ \item \code{none}: no pooling of K600 \item @@ -158,7 +158,7 @@ #' mm_name('sim', err_proc_acor=TRUE) #' mm_name('bayes', pool_K600='binned') mm_name <- function( - type=c('mle','bayes','bayes_2station','night','Kmodel','sim'), + type=c('mle','bayes','bayes_2s','night','Kmodel','sim'), #pool_GPP='none', pool_ER='none', pool_eoi='alldays', pool_epc='alldays', pool_epi='alldays', pool_K600=c('none', 'normal','normal_sdzero','normal_sdfixed', @@ -177,26 +177,26 @@ mm_name <- function( deficit_src=c('DO_mod','DO_obs','DO_obs_filter','NA'), engine=c('stan','nlm','lm','mean','loess','rnorm'), check_validity=TRUE) { - - # determine type. 'bayes_2station' is matched exactly, before match.arg's + + # determine type. 'bayes_2s' is matched exactly, before match.arg's # partial-prefix matching, because 'b' would otherwise be an ambiguous - # abbreviation between 'bayes' and 'bayes_2station' -- so unlike the other - # types, 'bayes_2station' must be spelled out in full (no abbreviations). + # abbreviation between 'bayes' and 'bayes_2s' -- so unlike the other + # types, 'bayes_2s' must be spelled out in full (no abbreviations). # match.arg's choices are narrowed to exclude it so that pre-existing # abbreviations like 'b' (-> 'bayes') and 'm' (-> 'mle') stay unambiguous. if(missing(type)) { type <- eval(formals(mm_name)$type)[1] - } else if(length(type) == 1 && identical(type, 'bayes_2station')) { - type <- 'bayes_2station' + } else if(length(type) == 1 && identical(type, 'bayes_2s')) { + type <- 'bayes_2s' } else { - type <- match.arg(type, choices=setdiff(eval(formals(mm_name)$type), 'bayes_2station')) + type <- match.arg(type, choices=setdiff(eval(formals(mm_name)$type), 'bayes_2s')) } - # bayes_2station has a single fixed model structure rather than being built + # bayes_2s has a single fixed model structure rather than being built # from combinations of pool_K600/err_*/ode_method/GPP_fun/ER_fun/ # deficit_src/engine, so skip the argument-combination machinery below and # return the one valid name directly - if(type == 'bayes_2station') { + if(type == 'bayes_2s') { mmname <- 'b2_np_oi_tr_plrckm.stan' check_validity <- if(!is.logical(check_validity)) stop("need check_validity to be a logical of length 1") else check_validity[1] if(isTRUE(check_validity)) mm_validate_name(mmname) @@ -209,7 +209,7 @@ mm_name <- function( relevant_args <- names(formals(mm_name)) %>% .[!(. %in% c('type','check_validity'))] } else { # only one argument allowed for Kmodel - relevant_args <- 'engine' + relevant_args <- 'engine' # directly specify all the rest pool_K600='complete' pool_all='complete' @@ -230,9 +230,9 @@ mm_name <- function( assign(ms, default_args[[ms]]) } } - - # check arguments and throw errors as needed. these checks define the names - # that are possible to create; will be supplemented by call to mm_valid_names + + # check arguments and throw errors as needed. these checks define the names + # that are possible to create; will be supplemented by call to mm_valid_names # to see if a specific arg combo is actually implemented if(type != 'Kmodel') { pool_K600 <- match.arg(pool_K600) @@ -254,7 +254,7 @@ mm_name <- function( engine <- match.arg(engine) if(!(engine %in% list(bayes='stan', mle='nlm', night='lm', Kmodel=c('mean','lm','loess'), sim='rnorm')[[type]])) stop("mismatch between type (",type,") and engine (",engine,")") - + # make the name mmname <- paste0( c(bayes='b', mle='m', night='n', Kmodel='K', sim='s')[[type]], '_', @@ -262,19 +262,19 @@ mm_name <- function( c(none_or_fitted='', sdzero='0', sdfixed='x')[[tryCatch(strsplit(pool_K600, '_')[[1]][[2]], error=function(e) 'none_or_fitted')]], c(none='np', partial='', complete='')[[pool_all]], '_', if(err_obs_iid) 'oi', if(err_proc_acor) 'pc', if(err_proc_iid) 'pi', if(err_proc_GPP) 'pp', '_', - c(Euler='Eu', pairmeans='pm', trapezoid='tr', rk2='r2', - lsoda='o1', lsode='o2', lsodes='o3', lsodar='o4', vode='o5', daspk='o6', euler='eu', rk4='o8', + c(Euler='Eu', pairmeans='pm', trapezoid='tr', rk2='r2', + lsoda='o1', lsode='o2', lsodes='o3', lsodar='o4', vode='o5', daspk='o6', euler='eu', rk4='o8', ode23='o9', ode45='o10', radau='o11', bdf='o12', bdf_d='o13', adams='o14', impAdams='o15', impAdams_d='o16', 'NA'='')[[ode_method]], '_', c(linlight='pl', satlight='ps', satlightq10temp='pq', 'NA'='')[[GPP_fun]], c(constant='rc', q10temp='rq', 'NA'='')[[ER_fun]], - c(DO_mod='km', DO_obs='ko', DO_obs_filter='kf', 'NA'='')[[deficit_src]], + c(DO_mod='km', DO_obs='ko', DO_obs_filter='kf', 'NA'='')[[deficit_src]], '.', engine) - + # check validity if requested check_validity <- if(!is.logical(check_validity)) stop("need check_validity to be a logical of length 1") else check_validity[1] if(isTRUE(check_validity)) mm_validate_name(mmname) - + # return mmname } diff --git a/R/mm_parse_name.R b/R/mm_parse_name.R index 1d8d824..d169369 100644 --- a/R/mm_parse_name.R +++ b/R/mm_parse_name.R @@ -1,20 +1,20 @@ #' Parse a model name into its features -#' -#' Returns a data.frame with one column per model structure detail and one row -#' per `model_name` supplied to this function. See \code{?\link{mm_name}} for a +#' +#' Returns a data.frame with one column per model structure detail and one row +#' per `model_name` supplied to this function. See \code{?\link{mm_name}} for a #' description of each of the data.frame columns that is returned. -#' -#' Custom model files (for MCMC) may have additional characters after an -#' underscore at the end of the name and before the prefix. For example, +#' +#' Custom model files (for MCMC) may have additional characters after an +#' underscore at the end of the name and before the prefix. For example, #' 'b_np_pcpi_eu_ko.stan' and 'b_np_pcpi_eu_ko_v2.stan' are parsed the same; the #' _v2 is ignored by this function. -#' +#' #' @seealso The converse of this function is \code{\link{mm_name}}. -#' +#' #' @param model_name character: the model name #' @param expand logical: should additional columns such as model_name and #' pool_K600_type be added? If expand=TRUE then the result cannot be passed -#' directly back into mm_name, but the additional columns may be helpful for +#' directly back into mm_name, but the additional columns may be helpful for #' interpreting the model structure. #' @import dplyr #' @importFrom stats na.omit @@ -25,25 +25,24 @@ mm_parse_name <- function(model_name, expand=FALSE) { # define function that gets used to parse prk_terms - match_or_NA <- function(key, pairs) { - matches <- c(unname(na.omit(key[pairs]))) + match_or_NA <- function(key, pairs) { + matches <- c(unname(na.omit(key[pairs]))) if(length(matches) == 0) { 'NA' } else if(length(matches) > 1) { - stop('found too many matches in PRK terms') + stop('found too many matches in PRK terms') } else { matches } } # parse the name parsed <- strsplit(basename(model_name), "_|\\.") sapply(1:length(parsed), function(pnum) if(length(parsed[[pnum]]) <= 5) stop('missing one or more pieces in name: ', model_name[pnum])) - # 'b2' (not just 'b') is the whole first token for two-station model files - # (e.g., 'b2_np_oi_tr_plrckm.stan'), since strsplit on "_|\\." above never - # splits within a token; an exact-match entry here is therefore sufficient - # and requires no change to the token-extraction logic itself - type <- unname(c(b='bayes', b2='bayes_2station', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) + # the "_|\\." split regex above handles 'b2' correctly. No change + # to the token-extraction logic was needed to support two-station names -- + # only a new lookup entry below, mapping the 'b2' token to its type name. + type <- unname(c(b='bayes', b2='bayes_2s', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) pool_K600 <- unname(c( - np='none', + np='none', Kn='normal', Kn0='normal_sdzero', Knx='normal_sdfixed', Kl='linear', Kl0='linear_sdzero', Klx='linear_sdfixed', Kb='binned', Kb0='binned_sdzero', Kbx='binned_sdfixed', @@ -63,8 +62,8 @@ mm_parse_name <- function(model_name, expand=FALSE) { err_proc_iid <- grepl('pi', sapply(parsed, `[`, 3)) err_proc_GPP <- grepl('pp', sapply(parsed, `[`, 3)) ode_method <- unname( - c(Eu='Euler', pm='pairmeans', tr='trapezoid', r2='rk2', o1='lsoda', o2='lsode', o3='lsodes', - o4='lsodar', o5='vode', o6='daspk', o7='euler', eu='euler', o8='rk4', o9='ode23', o10='ode45', o11='radau', + c(Eu='Euler', pm='pairmeans', tr='trapezoid', r2='rk2', o1='lsoda', o2='lsode', o3='lsodes', + o4='lsodar', o5='vode', o6='daspk', o7='euler', eu='euler', o8='rk4', o9='ode23', o10='ode45', o11='radau', o12='bdf', o13='bdf_d', o14='adams', o15='impAdams', o16='impAdams_d')[sapply(parsed, `[`, 4)]) prk_terms <- bind_rows(lapply(parsed, function(parsed1) { prk_term <- parsed1[5] @@ -79,7 +78,7 @@ mm_parse_name <- function(model_name, expand=FALSE) { ER_fun <- prk_terms$ER_fun deficit_src <- prk_terms$deficit_src engine <- sapply(parsed, function(vec) vec[length(vec)]) # the last one - leaves room for custom name endings before the suffix - + # combine the parsed pieces into a data.frame df <- data.frame( model_name=model_name, @@ -95,10 +94,10 @@ mm_parse_name <- function(model_name, expand=FALSE) { GPP_fun=ifelse(is.na(GPP_fun), 'NA', GPP_fun), ER_fun=ifelse(is.na(ER_fun), 'NA', ER_fun), deficit_src=ifelse(is.na(deficit_src), 'NA', deficit_src), - engine=ifelse(is.na(engine), 'NA', engine), + engine=ifelse(is.na(engine), 'NA', engine), stringsAsFactors=FALSE) - + if(!expand) df$model_name <- df$pool_K600_type <- df$pool_K600_sd <- NULL - + df } diff --git a/R/mm_ts_prep_data.R b/R/mm_ts_prep_data.R deleted file mode 100644 index 5414261..0000000 --- a/R/mm_ts_prep_data.R +++ /dev/null @@ -1,125 +0,0 @@ -#' Reshape long-format two-station data into the list expected by the -#' two-station Stan model -#' -#' Time-shifts the upstream DO series to match the travel time between -#' stations, then pivots the result into the \code{n_obs x n_days} matrices -#' expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} -#' (see \code{\link{metab_2station}}). -#' -#' The upstream observation that "matches" a given downstream observation at -#' row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / -#' timestep_days)} timesteps earlier, where \code{timestep_days} is the -#' median timestep of \code{data$solar.time}. This must be computed the same -#' way as in \code{\link{metab_2station}}'s lead-in check, so that the -#' \code{max(lag)} computed here always agrees with the lead-in requirement -#' already validated there. The first \code{max(lag)} rows of \code{data} are -#' lead-in rows: they supply upstream DO for the shift but are never -#' themselves treated as modeled (downstream) observations. -#' -#' @param data data.frame as validated by \code{\link{mm_validate_data}} for -#' \code{\link{metab_2station}}: must contain \code{solar.time}, -#' \code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, -#' \code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, -#' \code{travel.time}, sorted ascending by \code{solar.time}, and must -#' include the lead-in rows required to cover the longest travel time (see -#' \code{\link{metab_2station}}). -#' @param specs optional list of model specs. If it contains -#' \code{K600_lnorm_meanlog} and/or \code{K600_lnorm_sdlog}, those values -#' are used; otherwise placeholder defaults of \code{log(3.48)} and -#' \code{0.5} are used. -#' @return a named list with all variables in the Stan model's data block: -#' \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, -#' \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, -#' \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} -#' matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} -#' @importFrom unitted v -#' @export -mm_ts_prep_data <- function(data, specs=NULL) { - - # strip units; Stan cannot handle unitted vectors/matrices - data <- v(data) - - # timestep_days must match metab_2station()'s lead-in check exactly (median - # timestep, in days), so that max_lag here agrees with what was validated - # there - timestep_days <- stats::median(as.numeric(diff(data$solar.time), units='days')) - - # lag, in timesteps, that the upstream series must be shifted by to line up - # with each row's downstream observation - lag <- round(data$travel.time / timestep_days) - max_lag <- max(lag) - n_total <- nrow(data) - if(n_total <= max_lag) { - stop( - 'not enough lead-in rows to cover the longest travel.time (', max_lag, ' timesteps ', - 'implied, but only ', n_total, ' rows supplied); this should already have been caught ', - 'by metab_2station()') - } - - # the first max_lag rows are lead-in only (upstream DO used for the shift, - # but never modeled themselves). every row i > max_lag is guaranteed to - # have a valid shift target (i - lag[i] >= 1) because lag[i] <= max_lag - keep <- seq.int(max_lag + 1, n_total) - shift_idx <- keep - lag[keep] - if(any(shift_idx < 1)) { - # should be unreachable given max_lag's definition; guards against - # programming errors rather than expected user input - stop('internal error: upstream shift index falls before the first row of data') - } - - modeled <- data.frame( - solar.time = data$solar.time[keep], - DO_obs_up = data$DO.obs.up[shift_idx], - DO_sat_up = data$DO.sat.up[shift_idx], - DO_obs_down = data$DO.obs.down[keep], - DO_sat_down = data$DO.sat.down[keep], - light = data$light[keep], - depth = data$depth[keep], - temp_water = data$temp.water[keep], - travel_time = data$travel.time[keep] - ) - - # pivot into n_obs x n_days matrices, one column per unique date, following - # the same time_by_date_matrix approach used for the 1-station Stan models - # (see prepdata_bayes() in metab_bayes.R) - date_vec <- as.character(as.Date(modeled$solar.time)) - date_table <- table(date_vec) - n_days <- length(date_table) - n_obs_per_day <- unique(unname(date_table)) - if(length(n_obs_per_day) > 1) { - stop( - 'dates have differing numbers of modeled rows after lead-in removal; ', - 'observations cannot be combined into a matrix: ', - paste(sprintf('%s (%d rows)', names(date_table), date_table), collapse=', ')) - } - n_obs <- n_obs_per_day - - to_matrix <- function(vec) matrix(vec, nrow=n_obs, ncol=n_days, byrow=FALSE) - - # confirm each date occupies a contiguous block of rows, i.e., that data - # was sorted by solar.time; otherwise the matrix pivot below would silently - # scramble which rows belong to which date - date_mat <- to_matrix(date_vec) - unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) - if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { - stop('data must be sorted by solar.time so that each date occupies a contiguous block of rows') - } - - K600_lnorm_meanlog <- if(!is.null(specs$K600_lnorm_meanlog)) specs$K600_lnorm_meanlog else log(3.48) - K600_lnorm_sdlog <- if(!is.null(specs$K600_lnorm_sdlog)) specs$K600_lnorm_sdlog else 0.5 - - list( - n_obs = n_obs, - n_days = n_days, - DO_obs_up = to_matrix(modeled$DO_obs_up), - DO_sat_up = to_matrix(modeled$DO_sat_up), - DO_obs_down = to_matrix(modeled$DO_obs_down), - DO_sat_down = to_matrix(modeled$DO_sat_down), - light = to_matrix(modeled$light), - depth = to_matrix(modeled$depth), - temp_water = to_matrix(modeled$temp_water), - travel_time = to_matrix(modeled$travel_time), - K600_lnorm_meanlog = K600_lnorm_meanlog, - K600_lnorm_sdlog = K600_lnorm_sdlog - ) -} diff --git a/R/mm_valid_names.R b/R/mm_valid_names.R index 1a98d1b..9423984 100644 --- a/R/mm_valid_names.R +++ b/R/mm_valid_names.R @@ -10,7 +10,7 @@ #' @examples #' mm_valid_names('mle') #' @export -mm_valid_names <- function(type=c('bayes','bayes_2station','mle','night','Kmodel','sim')) { +mm_valid_names <- function(type=c('bayes','bayes_2s','mle','night','Kmodel','sim')) { type <- match.arg(type, several.ok=TRUE) @@ -39,7 +39,7 @@ mm_valid_names <- function(type=c('bayes','bayes_2station','mle','night','Kmodel mnames <- grep('^b_', dir(system.file('models', package='streamMetabolizer')), value=TRUE) favorites <- c('b_np_oipi_tr_plrckm.stan','b_np_oi_tr_plrckm.stan','b_np_pi_tr_plrckm.stan','b_np_oipp_tr_plrckm.stan') }, - bayes_2station={ + bayes_2s={ # single fixed model structure; no combinatorial name-building needed mnames <- 'b2_np_oi_tr_plrckm.stan' favorites <- mnames diff --git a/data-raw/two_station_example.R b/data-raw/two_station_example.R index 0142b89..b2749eb 100644 --- a/data-raw/two_station_example.R +++ b/data-raw/two_station_example.R @@ -1,10 +1,9 @@ -# Builds data/two_station_example.rda from the VFTS paper's published input -# data. Not run automatically as part of the package build/check; run +# Builds data/two_station_example.rda from the VFTS (Variable Flow +# Two-Station) paper's published input data. Not run automatically as part +# of the package build/check; run # manually (with the package root as the working directory) whenever the # example dataset needs to be regenerated. # -# File source: 2_station/Data/2_VFTS_and_One-station_model_input.csv, a sibling -# directory of the package root (not included in the package itself) # Download from ScienceBase: https://www.sciencebase.gov/catalog/item/6887d457d4be024722b4aae2 @@ -36,14 +35,14 @@ modeled_start <- as.POSIXct('2011-07-31 00:00:00', tz='UTC') modeled_end <- as.POSIXct('2011-08-29 23:45:00', tz='UTC') timestep_days <- 15/(24*60) -# metab_2station()'s upstream-DO lag shift (see mm_ts_prep_data() and -# metab_2station()'s "Two-station data requirements" section) needs +# metab_bayes_2s()'s upstream-DO lag shift (see prepdata_bayes_2s() and +# metab_bayes_2s()'s "Two-station data requirements" section) needs # max_lag = max(round(travel.time / timestep_days)) rows of lead-in -# immediately before modeled_start -- and because mm_ts_prep_data() trims +# immediately before modeled_start -- and because prepdata_bayes_2s() trims # max_lag rows off the *start of the whole array*, not off each calendar # day, that lead-in window must be exactly max_lag rows (not e.g. a whole # extra day) or the first modeled date ends up with a different row count -# than the rest, which mm_ts_prep_data() rejects. max_lag is computed from a +# than the rest, which prepdata_bayes_2s() rejects. max_lag is computed from a # generous 2-day candidate lead-in window and then trimmed to size. candidate_start <- modeled_start - as.difftime(2, units='days') candidate <- vfts2 %>% filter(datetime >= candidate_start, datetime <= modeled_end) @@ -54,7 +53,7 @@ vfts2_window <- vfts2 %>% filter(datetime >= lead_in_start, datetime <= modeled_end) # confirm the window is gap-free at the native 15-min timestep, and that -# trimming the lead-in rows (as mm_ts_prep_data() does) leaves exactly 30 +# trimming the lead-in rows (as prepdata_bayes_2s() does) leaves exactly 30 # modeled dates with equal row counts stopifnot(all(abs(diff(as.numeric(vfts2_window$datetime)) - 15*60) < 1e-6)) modeled_dates <- as.Date(vfts2_window$datetime[seq.int(max_lag+1, nrow(vfts2_window))]) diff --git a/man/get_params.Rd b/man/get_params.Rd index 644e237..2a5b3c8 100644 --- a/man/get_params.Rd +++ b/man/get_params.Rd @@ -1,10 +1,12 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_Kmodel.R, -% R/metab_bayes.R, R/metab_model.get_params.R, R/metab_sim.R +% R/metab_bayes.R, R/metab_bayes_2s.R, R/metab_model.get_params.R, +% R/metab_sim.R \name{get_params} \alias{get_params} \alias{get_params.metab_Kmodel} \alias{get_params.metab_bayes} +\alias{get_params.metab_bayes_2s} \alias{get_params.metab_model} \alias{get_params.metab_sim} \title{Extract the metabolism parameters (fitted and/or fixed) from a model.} @@ -42,6 +44,15 @@ get_params( attach.units = deprecated() ) +\method{get_params}{metab_bayes_2s}( + metab_model, + date_start = NA, + date_end = NA, + uncertainty = c("sd", "ci", "none"), + messages = TRUE, + ... +) + \method{get_params}{metab_model}( metab_model, date_start = NA, @@ -109,22 +120,34 @@ parameters describing the rates and/or shapes of GPP, ER, or reaeration. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_Kmodel}: Make daily re-predictions of K600.daily based on the +\item \code{get_params(metab_Kmodel)}: Make daily re-predictions of K600.daily based on the across-days model of K600.daily versus predictors. Only returns estimates for K600.daily, not any of the other daily parameters -\item \code{metab_bayes}: Does a little formatting to convert from Stan output +\item \code{get_params(metab_bayes)}: Does a little formatting to convert from Stan output to streamMetabolizer parameter names; otherwise the same as \code{get_params.metab_model} -\item \code{metab_model}: This implementation is shared by many model types - -\item \code{metab_sim}: Generates new simulated values for daily parameters if +\item \code{get_params(metab_bayes_2s)}: Does the same Stan-output-to-streamMetabolizer +renaming as \code{get_params.metab_bayes}, but (unlike that method) does +not delegate the rest of the work to \code{get_params.metab_model} via +\code{NextMethod()}: that generic implementation looks up parameter +names via \code{get_param_names()}, which builds an ODE-based dDOdt +function using streamMetabolizer's one-station instantaneous-rate +framework (\code{ode_method}/\code{GPP_fun}/\code{ER_fun}/ +\code{deficit_src}) -- machinery that doesn't apply to the two-station +steady-state model's daily GPP/ER/K600 parameters. \code{fixed} +column/star annotations (relevant only to models that can take fixed +daily parameters from \code{data_daily}) are not supported here. + +\item \code{get_params(metab_model)}: This implementation is shared by many model types + +\item \code{get_params(metab_sim)}: Generates new simulated values for daily parameters if they were described with evaluatable expressions in \code{\link{specs}}, or returns the fixed values for daily parameters if they were set in \code{data_daily} -}} +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -135,10 +158,10 @@ get_params(mm, date_start=get_fit(mm)$date[2]) \code{\link{predict_metab}} for daily average rates of GPP and ER Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_specs}()}, diff --git a/man/metab_Kmodel-class.Rd b/man/metab_Kmodel-class.Rd index 74fe446..4036c68 100644 --- a/man/metab_Kmodel-class.Rd +++ b/man/metab_Kmodel-class.Rd @@ -12,6 +12,7 @@ available data to reach better, less variable daily estimates of K \seealso{ Other metab.model.classes: \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, diff --git a/man/metab_Kmodel.Rd b/man/metab_Kmodel.Rd index f8368f6..6f2f37b 100644 --- a/man/metab_Kmodel.Rd +++ b/man/metab_Kmodel.Rd @@ -135,6 +135,7 @@ plot_metab_preds(mm3) \seealso{ Other metab_model: \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}}, \code{\link{metab_sim}} diff --git a/man/metab_bayes-class.Rd b/man/metab_bayes-class.Rd index f74b4d3..7435bf4 100644 --- a/man/metab_bayes-class.Rd +++ b/man/metab_bayes-class.Rd @@ -11,6 +11,7 @@ and K for a given DO curve. \seealso{ Other metab.model.classes: \code{\link{metab_Kmodel-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, diff --git a/man/metab_bayes.Rd b/man/metab_bayes.Rd index e772795..ae1b3d1 100644 --- a/man/metab_bayes.Rd +++ b/man/metab_bayes.Rd @@ -80,6 +80,7 @@ file.edit(get_specs(mm)$model_path) \seealso{ Other metab_model: \code{\link{metab_Kmodel}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}}, \code{\link{metab_sim}} diff --git a/man/metab_2station-class.Rd b/man/metab_bayes_2s-class.Rd similarity index 75% rename from man/metab_2station-class.Rd rename to man/metab_bayes_2s-class.Rd index ca14b20..9051296 100644 --- a/man/metab_2station-class.Rd +++ b/man/metab_bayes_2s-class.Rd @@ -1,11 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/metab_2station.R +% Please edit documentation in R/metab_bayes_2s.R \docType{class} -\name{metab_2station-class} -\alias{metab_2station-class} -\title{Metabolism model fitted by two-station (VFTS) Bayesian MCMC} +\name{metab_bayes_2s-class} +\alias{metab_bayes_2s-class} +\title{Metabolism model fitted by two-station (Variable Flow Two-Station) Bayesian +MCMC} \description{ -\code{metab_2station} models use Bayesian MCMC methods to fit values of +\code{metab_bayes_2s} models use Bayesian MCMC methods to fit values of GPP, ER, and K600 from paired upstream/downstream DO curves. This class inherits from \code{metab_bayes} (same \code{log}/\code{mcmc}/ \code{mcmc_data}/\code{compile_time} slots, and therefore the same diff --git a/man/metab_2station.Rd b/man/metab_bayes_2s.Rd similarity index 65% rename from man/metab_2station.Rd rename to man/metab_bayes_2s.Rd index e7ba103..3e7fe83 100644 --- a/man/metab_2station.Rd +++ b/man/metab_bayes_2s.Rd @@ -1,11 +1,11 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/metab_2station.R -\name{metab_2station} -\alias{metab_2station} -\title{Two-station Bayesian metabolism model fitting function (stub)} +% Please edit documentation in R/metab_bayes_2s.R +\name{metab_bayes_2s} +\alias{metab_bayes_2s} +\title{Two-station Bayesian metabolism model fitting function} \usage{ -metab_2station( - specs = specs(mm_name("bayes_2station")), +metab_bayes_2s( + specs = specs(mm_name("bayes_2s")), data = mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, depth, temp.water, travel.time), data_daily = mm_data(date, optional = "all"), @@ -35,16 +35,25 @@ description.} within the metab_model object} } \value{ -Not yet implemented; currently always errors after data validation. +A metab_bayes_2s object containing the fitted model. This object + can be inspected with the functions in the + \code{\link{metab_model_interface}} and also \code{\link{get_mcmc}}. } \description{ -Fits a two-station (upstream/downstream, a.k.a. VFTS) Bayesian model to -estimate GPP and ER from paired upstream and downstream DO, temperature, -light, and travel-time data. This function is currently a stub: it -validates the \code{data} argument and enforces two-station-specific data -requirements, but does not yet fit a model. See \code{\link{mm_name}} to -choose a Bayesian model and \code{\link{specs}} for relevant options for -the \code{specs} argument. +Fits a two-station (upstream/downstream, Variable Flow Two-Station) Bayesian +model to estimate GPP, ER, and K600 from paired upstream and downstream DO, +temperature, light, and travel-time data, using the single fixed Stan model +in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See \code{\link{mm_name}} to +choose a Bayesian model and \code{\link{specs}} for relevant options for the +\code{specs} argument. +} +\details{ +Unlike \code{\link{metab_bayes}}, which supports many model structures via +\code{split_dates}/\code{pool_K600}/etc., \code{metab_bayes_2s} always +fits every date jointly in a single Stan call (\code{specs$split_dates} is +forced to \code{FALSE} by \code{\link{specs}}), because the +upstream-downstream lag shift ties each date's first modeled rows to the +previous date's last rows. } \section{Two-station data requirements}{ In addition to the checks diff --git a/man/metab_mle-class.Rd b/man/metab_mle-class.Rd index 62f2301..4017dc1 100644 --- a/man/metab_mle-class.Rd +++ b/man/metab_mle-class.Rd @@ -12,6 +12,7 @@ likelihood to fit values of GPP, ER, and K for a given DO curve. Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_mle.Rd b/man/metab_mle.Rd index 1c9ab90..de2211d 100644 --- a/man/metab_mle.Rd +++ b/man/metab_mle.Rd @@ -72,6 +72,7 @@ plot_DO_preds(predict_DO(mm)) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_night}}, \code{\link{metab_sim}} } diff --git a/man/metab_model-class.Rd b/man/metab_model-class.Rd index e3070d9..7103041 100644 --- a/man/metab_model-class.Rd +++ b/man/metab_model-class.Rd @@ -34,6 +34,7 @@ function.} Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_night-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_night-class.Rd b/man/metab_night-class.Rd index dcf0660..362b14d 100644 --- a/man/metab_night-class.Rd +++ b/man/metab_night-class.Rd @@ -12,6 +12,7 @@ given DO time series. Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_night.Rd b/man/metab_night.Rd index e2fa158..8aa9db3 100644 --- a/man/metab_night.Rd +++ b/man/metab_night.Rd @@ -57,6 +57,7 @@ plot_DO_preds(predict_DO(mm)) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_sim}} } diff --git a/man/metab_sim-class.Rd b/man/metab_sim-class.Rd index 59c5922..eb74560 100644 --- a/man/metab_sim-class.Rd +++ b/man/metab_sim-class.Rd @@ -12,6 +12,7 @@ including GPP, ER, and K600 values Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}} diff --git a/man/metab_sim.Rd b/man/metab_sim.Rd index fa6e0bf..3dde3a6 100644 --- a/man/metab_sim.Rd +++ b/man/metab_sim.Rd @@ -108,6 +108,7 @@ library(ggplot2) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}} } diff --git a/man/mm_generate_mcmc_file.Rd b/man/mm_generate_mcmc_file.Rd index 144d924..1884b17 100644 --- a/man/mm_generate_mcmc_file.Rd +++ b/man/mm_generate_mcmc_file.Rd @@ -23,13 +23,16 @@ mm_generate_mcmc_file( \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} \item{pool_K600}{character. [How] should the model pool information among days to get more consistent daily estimates for K600? Options (see Details diff --git a/man/mm_name.Rd b/man/mm_name.Rd index abf4a4d..c1e8037 100644 --- a/man/mm_name.Rd +++ b/man/mm_name.Rd @@ -5,7 +5,7 @@ \title{Find the name of a model by its features} \usage{ mm_name( - type = c("mle", "bayes", "night", "Kmodel", "sim"), + type = c("mle", "bayes", "bayes_2s", "night", "Kmodel", "sim"), pool_K600 = c("none", "normal", "normal_sdzero", "normal_sdfixed", "linear", "linear_sdzero", "linear_sdfixed", "binned", "binned_sdzero", "binned_sdfixed", "complete"), @@ -27,13 +27,16 @@ mm_name( \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} \item{pool_K600}{character. [How] should the model pool information among days to get more consistent daily estimates for K600? Options (see Details diff --git a/man/mm_valid_names.Rd b/man/mm_valid_names.Rd index 0982a1e..fb578cd 100644 --- a/man/mm_valid_names.Rd +++ b/man/mm_valid_names.Rd @@ -4,19 +4,22 @@ \alias{mm_valid_names} \title{Get the valid names for a given model type or types} \usage{ -mm_valid_names(type = c("bayes", "mle", "night", "Kmodel", "sim")) +mm_valid_names(type = c("bayes", "bayes_2s", "mle", "night", "Kmodel", "sim")) } \arguments{ \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} } \description{ Returns a vector of the \code{model_name}s for the type[s] indicated. If diff --git a/man/predict_DO.Rd b/man/predict_DO.Rd index e7ba5e8..db2c026 100644 --- a/man/predict_DO.Rd +++ b/man/predict_DO.Rd @@ -1,9 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_Kmodel.R, -% R/metab_model.predict_DO.R, R/metab_night.R, R/metab_sim.R +% R/metab_bayes_2s.R, R/metab_model.predict_DO.R, R/metab_night.R, +% R/metab_sim.R \name{predict_DO} \alias{predict_DO} \alias{predict_DO.metab_Kmodel} +\alias{predict_DO.metab_bayes_2s} \alias{predict_DO.metab_model} \alias{predict_DO.metab_night} \alias{predict_DO.metab_sim} @@ -20,6 +22,8 @@ predict_DO( \method{predict_DO}{metab_Kmodel}(metab_model, date_start = NA, date_end = NA, ..., use_saved = TRUE) +\method{predict_DO}{metab_bayes_2s}(metab_model, date_start = NA, date_end = NA, ..., use_saved = TRUE) + \method{predict_DO}{metab_model}( metab_model, date_start = NA, @@ -64,22 +68,32 @@ oxygen. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_Kmodel}: Throws an error because models of type 'Kmodel' can't +\item \code{predict_DO(metab_Kmodel)}: Throws an error because models of type 'Kmodel' can't predict DO. \code{metab_Kmodel} predicts K at daily timesteps and usually knows nothing about GPP or ER. So it's not possible to predict DO from this model. Try passing the output to metab_mle and THEN predicting DO. -\item \code{metab_model}: This implementation is shared by many model types +\item \code{predict_DO(metab_bayes_2s)}: Two-station (Variable Flow Two-Station) models. +Returns a data.frame with columns \code{solar.time}, \code{DO.obs.down} +(the observed downstream DO from the input data), and \code{DO.mod.down} +(the posterior median of the two-station Stan model's fitted downstream DO) +-- unlike the one-station \code{predict_DO} methods, which return +\code{DO.obs}/\code{DO.mod}. The values are those computed once at fitting +time (see \code{\link{metab_bayes_2s}}); \code{use_saved=FALSE} (on-demand +recomputation from the fitted daily GPP/ER/K600 medians) is not +implemented. + +\item \code{predict_DO(metab_model)}: This implementation is shared by many model types -\item \code{metab_night}: Generate nighttime dissolved oxygen predictions from a +\item \code{predict_DO(metab_night)}: Generate nighttime dissolved oxygen predictions from a nighttime regression model. \code{metab_night} only fits ER and K, and only for the darkness hours, so predictions are only generated for those hours. -\item \code{metab_sim}: Simulate values for DO.obs (with process and +\item \code{predict_DO(metab_sim)}: Simulate values for DO.obs (with process and observation error), DO.mod (with process error only), and DO.pure (with no error). The errors are randomly generated on every new call to predict_DO. -}} +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -88,10 +102,10 @@ head(preds) } \seealso{ Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_params}()}, diff --git a/man/predict_metab.Rd b/man/predict_metab.Rd index ca653ff..f53bf06 100644 --- a/man/predict_metab.Rd +++ b/man/predict_metab.Rd @@ -1,9 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_bayes.R, -% R/metab_model.predict_metab.R +% R/metab_bayes_2s.R, R/metab_model.predict_metab.R \name{predict_metab} \alias{predict_metab} \alias{predict_metab.metab_bayes} +\alias{predict_metab.metab_bayes_2s} \alias{predict_metab.metab_model} \title{Predict metabolism from a fitted model.} \usage{ @@ -26,6 +27,8 @@ predict_metab( attach.units = deprecated() ) +\method{predict_metab}{metab_bayes_2s}(metab_model, date_start = NA, date_end = NA, ...) + \method{predict_metab}{metab_model}( metab_model, date_start = NA, @@ -91,16 +94,19 @@ GPP, ER, and K600. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_bayes}: Pulls daily metabolism estimates out of the Stan +\item \code{predict_metab(metab_bayes)}: Pulls daily metabolism estimates out of the Stan model results; looks for \code{GPP} or \code{GPP_daily} and for \code{ER} or \code{ER_daily} among the \code{params_out} (see \code{\link{specs}}), which means you can save just one (or both) of those sets of daily parameters when running the Stan model. Saving fewer parameters can help models run faster and use less RAM. -\item \code{metab_model}: This implementation is shared by many model types -}} +\item \code{predict_metab(metab_bayes_2s)}: Pulls daily GPP, ER, and K600 estimates out of +the two-station Stan model results. +\item \code{predict_metab(metab_model)}: This implementation is shared by many model types + +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -109,10 +115,10 @@ predict_metab(mm, date_start=get_fit(mm)$date[2]) } \seealso{ Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_params}()}, diff --git a/man/mm_ts_prep_data.Rd b/man/prepdata_bayes_2s.Rd similarity index 67% rename from man/mm_ts_prep_data.Rd rename to man/prepdata_bayes_2s.Rd index c85f226..8555719 100644 --- a/man/mm_ts_prep_data.Rd +++ b/man/prepdata_bayes_2s.Rd @@ -1,25 +1,28 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mm_ts_prep_data.R -\name{mm_ts_prep_data} -\alias{mm_ts_prep_data} +% Please edit documentation in R/metab_bayes_2s.R +\name{prepdata_bayes_2s} +\alias{prepdata_bayes_2s} \title{Reshape long-format two-station data into the list expected by the two-station Stan model} \usage{ -mm_ts_prep_data(data, specs = NULL) +prepdata_bayes_2s(data, specs = NULL) } \arguments{ \item{data}{data.frame as validated by \code{\link{mm_validate_data}} for -\code{\link{metab_2station}}: must contain \code{solar.time}, +\code{\link{metab_bayes_2s}}: must contain \code{solar.time}, \code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, \code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, \code{travel.time}, sorted ascending by \code{solar.time}, and must include the lead-in rows required to cover the longest travel time (see -\code{\link{metab_2station}}).} +\code{\link{metab_bayes_2s}}).} -\item{specs}{optional list of model specs. If it contains -\code{K600_lnorm_meanlog} and/or \code{K600_lnorm_sdlog}, those values -are used; otherwise placeholder defaults of \code{log(3.48)} and -\code{0.5} are used.} +\item{specs}{a list of model specs (see \code{\link{specs}}), expected to +already contain \code{K600_lnorm_meanlog} and \code{K600_lnorm_sdlog} +-- e.g., the object returned by \code{specs(mm_name('bayes_2s'))}, which +populates them with sensible defaults. This function does not supply +its own fallback values; if \code{specs} is omitted or missing these +fields, the resulting Stan data will contain NULL/missing values for +them.} } \value{ a named list with all variables in the Stan model's data block: @@ -32,16 +35,17 @@ a named list with all variables in the Stan model's data block: Time-shifts the upstream DO series to match the travel time between stations, then pivots the result into the \code{n_obs x n_days} matrices expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} -(see \code{\link{metab_2station}}). +(see \code{\link{metab_bayes_2s}}). } \details{ The upstream observation that "matches" a given downstream observation at row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / timestep_days)} timesteps earlier, where \code{timestep_days} is the median timestep of \code{data$solar.time}. This must be computed the same -way as in \code{\link{metab_2station}}'s lead-in check, so that the +way as in \code{\link{mm_validate_data_2station}}'s lead-in check, so that the \code{max(lag)} computed here always agrees with the lead-in requirement already validated there. The first \code{max(lag)} rows of \code{data} are lead-in rows: they supply upstream DO for the shift but are never themselves treated as modeled (downstream) observations. } +\keyword{internal} diff --git a/man/two_station_example.Rd b/man/two_station_example.Rd index 8bc0712..da4ab53 100644 --- a/man/two_station_example.Rd +++ b/man/two_station_example.Rd @@ -6,7 +6,7 @@ \title{Example two-station (VFTS) input data} \format{ A data.frame with 2904 rows and the 9 columns expected by - \code{\link{metab_2station}}'s \code{data} argument, each carrying + \code{\link{metab_bayes_2s}}'s \code{data} argument, each carrying \code{\link[unitted]{unitted}} units matching \code{\link{mm_data}}: \describe{ \item{solar.time}{POSIXct timestamp, UTC} @@ -43,13 +43,13 @@ two_station_example } \description{ A 30-day example dataset for fitting a two-station (upstream/downstream, -a.k.a. VFTS) metabolism model with \code{\link{metab_2station}}. It is a -subset of the \code{VFTS-2} (variable-travel-time) run from the published -two-station metabolism modeling dataset for a reach of the Colorado River -in Glen Canyon, covering 2011-07-31 through 2011-08-29 at the source -data's native 15-minute timestep, plus a lead-in block of upstream DO -observations (exactly as long as the longest travel time in the dataset -requires -- see \code{\link{metab_2station}}'s "Two-station data -requirements" section) immediately before 2011-07-31. +Variable Flow Two-Station) metabolism model with +\code{\link{metab_bayes_2s}}. It is a subset of the \code{VFTS-2} +(variable-travel-time) run from the published two-station metabolism modeling +dataset for a reach of the Colorado River in Glen Canyon, covering 2011-07-31 +through 2011-08-29 at the source data's native 15-minute timestep, plus a +lead-in block of upstream DO observations (exactly as long as the longest +travel time in the dataset requires -- see \code{\link{metab_bayes_2s}}'s +"Two-station data requirements" section) immediately before 2011-07-31. } \keyword{datasets} diff --git a/tests/testthat/test-metab_2station.R b/tests/testthat/test-metab_bayes_2s.R similarity index 82% rename from tests/testthat/test-metab_2station.R rename to tests/testthat/test-metab_bayes_2s.R index 18765c3..fdb541b 100644 --- a/tests/testthat/test-metab_2station.R +++ b/tests/testthat/test-metab_bayes_2s.R @@ -19,33 +19,33 @@ make_2station_data <- function(n=10, timestep_min=5, travel_time=0.01) { test_that("mm_validate_data catches missing required columns", { dat <- dplyr::select(make_2station_data(), -DO.obs.up) - expect_error(metab_2station(data=dat), "missing these columns") + expect_error(metab_bayes_2s(data=dat), "missing these columns") }) test_that("travel.time <= 0 triggers an error", { dat <- make_2station_data(travel_time=0) - expect_error(metab_2station(data=dat), "travel.time must be > 0") + expect_error(metab_bayes_2s(data=dat), "travel.time must be > 0") dat <- make_2station_data(travel_time=-0.01) - expect_error(metab_2station(data=dat), "travel.time must be > 0") + expect_error(metab_bayes_2s(data=dat), "travel.time must be > 0") }) test_that("travel.time >= 1 triggers an error with a units hint", { dat <- make_2station_data(travel_time=1) - expect_error(metab_2station(data=dat), "travel.time must be < 1 day.*incorrect units") + expect_error(metab_bayes_2s(data=dat), "travel.time must be < 1 day.*incorrect units") dat <- make_2station_data(travel_time=1.5) - expect_error(metab_2station(data=dat), "travel.time must be < 1 day.*incorrect units") + expect_error(metab_bayes_2s(data=dat), "travel.time must be < 1 day.*incorrect units") }) test_that("insufficient lead-in data triggers an error", { # 2 rows but max_lag=3 timesteps of upstream lead-in are needed dat <- make_2station_data(n=2) - expect_error(metab_2station(data=dat), "insufficient lead-in data") + expect_error(metab_bayes_2s(data=dat), "insufficient lead-in data") }) -# mm_ts_prep_data() ----------------------------------------------------- +# prepdata_bayes_2s() ----------------------------------------------------- # Build a two-day, unit-labeled data.frame with a known, traceable # DO.obs.up/DO.sat.up series (sequential integers) so the shift can be @@ -80,7 +80,7 @@ make_ts_data <- function(n_leadin=3, n_day1=10, n_day2=7, travel_time=0.01, unit test_that("upstream DO is shifted by the correct lag", { dat <- make_ts_data() - out <- mm_ts_prep_data(dat) + out <- prepdata_bayes_2s(dat) # max_lag=3, so modeled row i (original index i) uses upstream data from # original row (i - 3). day 1's 7 modeled rows are original rows 4:10, so @@ -94,7 +94,7 @@ test_that("upstream DO is shifted by the correct lag", { test_that("lead-in rows are excluded from the output matrices", { dat <- make_ts_data(n_leadin=3, n_day1=10, n_day2=7) - out <- mm_ts_prep_data(dat) + out <- prepdata_bayes_2s(dat) # 17 total rows in, 3 are lead-in-only, so 14 modeled rows should remain expect_equal(out$n_obs * out$n_days, nrow(dat) - 3) @@ -106,7 +106,7 @@ test_that("lead-in rows are excluded from the output matrices", { test_that("output matrices have n_obs x n_days dimensions", { dat <- make_ts_data() - out <- mm_ts_prep_data(dat) + out <- prepdata_bayes_2s(dat) expect_equal(out$n_obs, 7) expect_equal(out$n_days, 2) @@ -117,23 +117,25 @@ test_that("output matrices have n_obs x n_days dimensions", { test_that("all required Stan data block variables are present", { dat <- make_ts_data() - out <- mm_ts_prep_data(dat) + # K600_lnorm_meanlog/sdlog are owned by specs() (see PR D-6/I1); pass + # distinctive marker values here to confirm prepdata_bayes_2s() just reads + # them through from specs rather than computing its own defaults + out <- prepdata_bayes_2s(dat, specs=list(K600_lnorm_meanlog=1.23, K600_lnorm_sdlog=4.56)) expected_names <- c( 'n_obs','n_days','DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down', 'light','depth','temp_water','travel_time','K600_lnorm_meanlog','K600_lnorm_sdlog') expect_true(all(expected_names %in% names(out))) - # placeholder K600 lognormal priors, per PR D-3 - expect_equal(out$K600_lnorm_meanlog, log(3.48)) - expect_equal(out$K600_lnorm_sdlog, 0.5) + expect_equal(out$K600_lnorm_meanlog, 1.23) + expect_equal(out$K600_lnorm_sdlog, 4.56) }) test_that("units are stripped from all numeric outputs", { dat <- make_ts_data(unitted=TRUE) expect_true(is.unitted(dat)) - out <- mm_ts_prep_data(dat) + out <- prepdata_bayes_2s(dat, specs=list(K600_lnorm_meanlog=2.484907, K600_lnorm_sdlog=1.0)) for(varname in c('DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down','light','depth','temp_water','travel_time')) { expect_false(is.unitted(out[[varname]]), info=varname) } @@ -149,7 +151,7 @@ test_that("units are stripped from all numeric outputs", { test_that("mm_parse_name recognizes the b2_ prefix for two-station models", { parsed <- mm_parse_name('b2_np_oi_tr_plrckm.stan') - expect_equal(parsed$type, 'bayes_2station') + expect_equal(parsed$type, 'bayes_2s') # the rest of the name is shared syntax with one-station bayes models and # should parse the same way regardless of the b vs. b2 prefix expect_equal(parsed$pool_K600, 'none') @@ -168,18 +170,18 @@ test_that("mm_parse_name recognizes the b2_ prefix for two-station models", { }) -# mm_name() / mm_valid_names() / specs() for bayes_2station --------------- +# mm_name() / mm_valid_names() / specs() for bayes_2s --------------- -test_that("mm_name(type='bayes_2station') returns the single two-station model name", { - expect_equal(mm_name(type='bayes_2station'), 'b2_np_oi_tr_plrckm.stan') +test_that("mm_name(type='bayes_2s') returns the single two-station model name", { + expect_equal(mm_name(type='bayes_2s'), 'b2_np_oi_tr_plrckm.stan') }) -test_that("mm_valid_names('bayes_2station') returns the single two-station model name", { - expect_equal(mm_valid_names('bayes_2station'), 'b2_np_oi_tr_plrckm.stan') +test_that("mm_valid_names('bayes_2s') returns the single two-station model name", { + expect_equal(mm_valid_names('bayes_2s'), 'b2_np_oi_tr_plrckm.stan') }) -test_that("specs(mm_name('bayes_2station')) has the expected params_in/params_out/split_dates", { - sp <- specs(mm_name('bayes_2station')) +test_that("specs(mm_name('bayes_2s')) has the expected params_in/params_out/split_dates", { + sp <- specs(mm_name('bayes_2s')) expect_equal( sp$params_in, @@ -191,13 +193,13 @@ test_that("specs(mm_name('bayes_2station')) has the expected params_in/params_ou }) -# metab_2station() fitting, predict_metab(), predict_DO() ------------------ +# metab_bayes_2s() fitting, predict_metab(), predict_DO() ------------------ # Subset two_station_example to just a few modeled days for a faster test # fit. Naively slicing rows doesn't work: max_lag (the number of upstream # lead-in rows required) is recomputed from whatever travel.time values are # present in the slice, so an arbitrary row range can leave a partial first -# date once mm_ts_prep_data() trims max_lag rows off the front -- the same +# date once prepdata_bayes_2s() trims max_lag rows off the front -- the same # lead-in-sizing logic used in data-raw/two_station_example.R is needed here # too. subset_2station_data <- function(full_data, n_modeled_days) { @@ -224,11 +226,11 @@ test_that("metab() fits a two-station model and predict_metab()/predict_DO() wor small_dat <- subset_2station_data(two_station_example, n_modeled_days=3) sp <- specs( - mm_name('bayes_2station'), + mm_name('bayes_2s'), n_chains=1, n_cores=1, burnin_steps=100, saved_steps=100, verbose=FALSE) mm <- metab(specs=sp, data=small_dat) - expect_s4_class(mm, 'metab_2station') + expect_s4_class(mm, 'metab_bayes_2s') pm <- predict_metab(mm) expect_s3_class(pm, 'data.frame') From 2074de0e8a6f3f709a278ace772b5a635488e5d9 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 27 Jul 2026 13:58:30 -0700 Subject: [PATCH 09/17] Consolidate two-station validation into mm_validate_data_2station() 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. --- R/metab_bayes_2s.R | 13 +++--- R/mm_validate_data.R | 64 ++++++++++++++++++++++++++-- man/metab_bayes_2s.Rd | 14 +++--- man/mm_validate_data_2station.Rd | 22 ++++++++++ tests/testthat/test-metab_bayes_2s.R | 6 +-- 5 files changed, 103 insertions(+), 16 deletions(-) create mode 100644 man/mm_validate_data_2station.Rd diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 8b192d8..94461ed 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -30,11 +30,14 @@ utils::globalVariables(c(".", "metab_50pct", "DO.mod.down")) #' @section Two-station data requirements: In addition to the checks #' performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the #' reach travel time between the upstream and downstream stations, in days) -#' must be strictly positive and less than 1 day (values >= 1 usually -#' indicate travel time was supplied in the wrong units). There must also be -#' enough lead-in observations of upstream DO before the first row of -#' \code{data} to cover the longest travel time in the dataset, given the -#' (median) timestep of \code{data$solar.time}. +#' must be strictly positive and no greater than 8/24 days (8 hours). +#' Values above this limit either indicate travel time was supplied in the +#' wrong units (e.g., minutes or hours instead of days), or reflect a reach +#' whose actual travel time exceeds the 8-hour limit required to prevent +#' the previous day's light conditions from influencing the following +#' day's metabolism estimate. There must also be enough lead-in +#' observations of upstream DO before the first row of \code{data} to +#' cover the longest travel time in the dataset. #' #' @export #' @family metab_model diff --git a/R/mm_validate_data.R b/R/mm_validate_data.R index 401b9bc..0b7bf90 100644 --- a/R/mm_validate_data.R +++ b/R/mm_validate_data.R @@ -87,16 +87,74 @@ mm_validate_data <- function( data.units <- get_units(dat)[mismatched.units] expected.units <- get_units(expected.data)[mismatched.units] stop(paste0("unexpected units in ", data_type, ": ", paste0( - "(", 1:length(mismatched.units), ") ", + "(", 1:length(mismatched.units), ") ", names(data.units), " = ", data.units, ", expected ", expected.units, collapse="; ")), call.=FALSE) } } - + + # check travel.time bounds, if present. only two-station models supply + # this column, so this is a no-op for other model types. travel.time is + # expected in days; values outside (0, 8/24] either reflect a units + # mistake (e.g., minutes or hours rather than days) or a reach whose + # travel time exceeds the 8-hour limit needed to keep the previous + # day's light from bleeding into the following day's metabolism estimate + if('travel.time' %in% names(dat)) { + travel.time <- v(dat$travel.time) + if(any(travel.time <= 0)) { + stop('travel.time must be > 0', call.=FALSE) + } + if(any(travel.time > 8/24)) { + stop('travel.time must be <= 8/24 days (8 hours); values above this either suggest incorrect units ', + '(expected days, e.g. not minutes or hours) or a reach travel time that exceeds the 8-hour limit ', + "required to prevent the previous day's light conditions from influencing the following day's ", + 'metabolism estimate', call.=FALSE) + } + } + # return the data, whose columns may be reordered/filtered dat }) - + # return the data.frames, which may have had their columns reordered during validation and are packaged as a list return(dat_all) } + + +#' Two-station-specific data validation +#' +#' Checks the lead-in coverage requirement specific to +#' \code{\link{metab_bayes_2s}}: there must be enough lead-in observations of +#' upstream DO before the first modeled row to cover the longest travel time +#' in the dataset, given the (median) timestep of \code{data$solar.time}. +#' Column presence, timestamp validity, and travel.time bounds are expected +#' to have already been checked by \code{\link{mm_validate_data}}. +#' +#' @param data data.frame as returned by \code{\link{mm_validate_data}} for +#' \code{\link{metab_bayes_2s}}: must contain \code{solar.time} and +#' \code{travel.time}, sorted ascending by \code{solar.time}. +#' @keywords internal +mm_validate_data_2station <- function(data) { + + data_v <- v(data) + travel_time <- data_v$travel.time + solar_time <- data_v$solar.time + + # there must be enough lead-in rows of upstream DO before the first + # modeled row to cover the longest travel time in the dataset. timestep_days + # is the median observation interval, in days; max_lag is the number of + # timesteps by which upstream data must lead downstream predictions. The + # first max_lag rows of data serve only as lead-in and cannot themselves be + # modeled, so at least max_lag + 1 rows are required overall. + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + max_lag <- max(round(travel_time / timestep_days)) + if(nrow(data) <= max_lag) { + lead_in_needed <- max_lag - nrow(data) + 1 + stop(paste0( + 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', + max_lag, ' timestep(s), but only ', nrow(data), ' row(s) were supplied; ', + 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) + } + + invisible(NULL) +} diff --git a/man/metab_bayes_2s.Rd b/man/metab_bayes_2s.Rd index 3e7fe83..4f7571e 100644 --- a/man/metab_bayes_2s.Rd +++ b/man/metab_bayes_2s.Rd @@ -59,11 +59,15 @@ previous date's last rows. In addition to the checks performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the reach travel time between the upstream and downstream stations, in days) - must be strictly positive and less than 1 day (values >= 1 usually - indicate travel time was supplied in the wrong units). There must also be - enough lead-in observations of upstream DO before the first row of - \code{data} to cover the longest travel time in the dataset, given the - (median) timestep of \code{data$solar.time}. + must be strictly positive and no greater than 8/24 days (8 hours). + Values above this limit either indicate travel time was supplied in the + wrong units (e.g., minutes or hours instead of days), or reflect a reach + whose actual travel time exceeds the 8-hour limit required to prevent + the previous day's light conditions from influencing the following + day's metabolism estimate. There must also be enough lead-in + observations of upstream DO before the first row of \code{data} to + cover the longest travel time in the dataset, given the (median) + timestep of \code{data$solar.time}. } \seealso{ diff --git a/man/mm_validate_data_2station.Rd b/man/mm_validate_data_2station.Rd new file mode 100644 index 0000000..61984e3 --- /dev/null +++ b/man/mm_validate_data_2station.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_validate_data.R +\name{mm_validate_data_2station} +\alias{mm_validate_data_2station} +\title{Two-station-specific data validation} +\usage{ +mm_validate_data_2station(data) +} +\arguments{ +\item{data}{data.frame as returned by \code{\link{mm_validate_data}} for +\code{\link{metab_bayes_2s}}: must contain \code{solar.time} and +\code{travel.time}, sorted ascending by \code{solar.time}.} +} +\description{ +Checks the lead-in coverage requirement specific to +\code{\link{metab_bayes_2s}}: there must be enough lead-in observations of +upstream DO before the first modeled row to cover the longest travel time +in the dataset, given the (median) timestep of \code{data$solar.time}. +Column presence, timestamp validity, and travel.time bounds are expected +to have already been checked by \code{\link{mm_validate_data}}. +} +\keyword{internal} diff --git a/tests/testthat/test-metab_bayes_2s.R b/tests/testthat/test-metab_bayes_2s.R index fdb541b..1c8be70 100644 --- a/tests/testthat/test-metab_bayes_2s.R +++ b/tests/testthat/test-metab_bayes_2s.R @@ -30,12 +30,12 @@ test_that("travel.time <= 0 triggers an error", { expect_error(metab_bayes_2s(data=dat), "travel.time must be > 0") }) -test_that("travel.time >= 1 triggers an error with a units hint", { +test_that("travel.time > 8/24 days (8 hours) triggers an error with a units/limit hint", { dat <- make_2station_data(travel_time=1) - expect_error(metab_bayes_2s(data=dat), "travel.time must be < 1 day.*incorrect units") + expect_error(metab_bayes_2s(data=dat), "travel.time must be <= 8/24 days.*incorrect units") dat <- make_2station_data(travel_time=1.5) - expect_error(metab_bayes_2s(data=dat), "travel.time must be < 1 day.*incorrect units") + expect_error(metab_bayes_2s(data=dat), "travel.time must be <= 8/24 days.*incorrect units") }) test_that("insufficient lead-in data triggers an error", { From 8841358193d66c74d203560e2b79cb8280814f8f Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 27 Jul 2026 15:09:09 -0700 Subject: [PATCH 10/17] Move K600_lnorm_meanlog/sdlog defaults into specs() Matches the convention used by all other default spec values, rather than computing fallback defaults inline in prepdata_bayes_2s(). --- R/specs.R | 162 ++++++++++++++++++++++++++------------------------- man/specs.Rd | 24 +++++--- 2 files changed, 99 insertions(+), 87 deletions(-) diff --git a/R/specs.R b/R/specs.R index a1facb8..5528a63 100644 --- a/R/specs.R +++ b/R/specs.R @@ -250,9 +250,9 @@ #' proportional to light (with noise) and is applied to GPP rather than to #' dDO/dt. #' -#' @param K600_lnorm_meanlog hyperparameter for \code{type='bayes_2station'}. +#' @param K600_lnorm_meanlog hyperparameter for \code{type='bayes_2s'}. #' The mean of a lognormal prior distribution for K600_daily. -#' @param K600_lnorm_sdlog hyperparameter for \code{type='bayes_2station'}. The +#' @param K600_lnorm_sdlog hyperparameter for \code{type='bayes_2s'}. The #' standard deviation parameter of a lognormal prior distribution for #' K600_daily. #' @@ -343,39 +343,39 @@ #' specs(mm_name(type='bayes', pool_K600='normal')) #' @export specs <- function( - + ## All or several models - + model_name = mm_name(), engine, - + # inheritParams mm_model_by_ply day_start = 4, day_end = 28, - + # inheritParams mm_is_valid_day day_tests=c('full_day', 'even_timesteps', 'complete_data', 'pos_discharge', 'pos_depth'), required_timestep=NA, - - + + ## MLE - + # initial values - init.GPP.daily = 8, + init.GPP.daily = 8, init.Pmax = 10, init.alpha = 0.0001, - init.ER.daily = -10, + init.ER.daily = -10, init.ER20 = -10, init.K600.daily = 10, - - + + ## Bayes - + # model setup split_dates, keep_mcmcs = TRUE, keep_mcmc_data = TRUE, - + # hyperparameters for non-hierarchical GPP & ER GPP_daily_mu = 3.1, GPP_daily_lower = -Inf, @@ -387,41 +387,41 @@ specs <- function( ER_daily_mu = -7.1, ER_daily_upper = Inf, ER_daily_sigma = 7.1, - + # hyperparameters for non-hierarchical K600 K600_daily_meanlog = log(12), - + # hyperparameters for hierarchical K600 - normal K600_daily_meanlog_meanlog = log(12), K600_daily_meanlog_sdlog = 1.32, - + # hyperparameters for hierarchical K600 - linear. defaults should be # reasonably constrained, not too wide lnK600_lnQ_intercept_mu = 2, lnK600_lnQ_intercept_sigma = 2.4, lnK600_lnQ_slope_mu = 0, lnK600_lnQ_slope_sigma = 0.5, - - # hyperparameters for hierarchical K600 - binned. K600_daily ~ - # lognormal(K600_daily_nodes_meanlog[lnQ_bin], + + # hyperparameters for hierarchical K600 - binned. K600_daily ~ + # lognormal(K600_daily_nodes_meanlog[lnQ_bin], # K600_daily_nodes_sdlog[lnQ_bin]) with linear interpolation among bins before - # exponentiating. nodes_meanlog and nodes_sdlog may be length b = - # length(K600_daily_lnQ_nodes) or length 1 (to be replicated to length b). - # -8:6 covers almost all points in Raymond et al. 2012 and will therefore + # exponentiating. nodes_meanlog and nodes_sdlog may be length b = + # length(K600_daily_lnQ_nodes) or length 1 (to be replicated to length b). + # -8:6 covers almost all points in Raymond et al. 2012 and will therefore # always be too broad a range for a single stream. -3:3 will catch some # streams to rivers as a first cut, though users should still modify K600_lnQ_nodes_centers = -3:3, # the x=lnQ values for the nodes K600_lnQ_nodediffs_sdlog = 0.5, # for centers 1 apart; for centers 0.2 apart, use 1/5 of this K600_lnQ_nodes_meanlog = rep(log(12), length(K600_lnQ_nodes_centers)), # distribs for the y=K600 values of the nodes K600_lnQ_nodes_sdlog = rep(1.32, length(K600_lnQ_nodes_centers)), - + # hyperparameters for any K pooling or non-pooling strategy K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none=1, normal_sdfixed=0.05, NA), K600_daily_sigma = switch(mm_parse_name(model_name)$pool_K600, linear_sdfixed=10, binned_sdfixed=5, NA), K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal=0.05, NA), K600_daily_sigma_sigma = switch(mm_parse_name(model_name)$pool_K600, linear=1.2, binned=0.24, NA), # normal_sdzero, linear_sdzero, and binned_sdzero all have no parameters for this - + # hyperparameters for error terms err_obs_iid_sigma_scale = 0.03, err_proc_iid_sigma_scale = 5, @@ -430,14 +430,16 @@ specs <- function( err_proc_acor_sigma_scale = 1, err_mult_GPP_sdlog_sigma = 1, - # hyperparameters for two-station (bayes_2station) K600. GPP_daily_mu, - # GPP_daily_sigma, ER_daily_mu, and ER_daily_sigma above are reused as-is + # hyperparameters for two-station (bayes_2s) K600. GPP_daily_mu, + # GPP_daily_sigma, ER_daily_mu, and ER_daily_sigma above are reused as-is. + # These are the sole source of the K600 lognormal prior defaults -- + # prepdata_bayes_2s() reads them from specs$K600_lnorm_meanlog/sdlog K600_lnorm_meanlog = 2.484907, K600_lnorm_sdlog = 1.0, # vector of hyperparameters to include as MCMC data params_in, - + # inheritParams runstan_bayes params_out, n_chains = 4, @@ -446,22 +448,22 @@ specs <- function( saved_steps = 500, thin_steps = 1, verbose = FALSE, - - + + ## Kmodel - + #inheritParams prepdata_Kmodel weights = c("K600/CI"), # 'K600/CI' is argued for in stream_metab_usa issue #64 filters = c(CI.max=NA, discharge.daily.max=NA, velocity.daily.max=NA), - + #inheritParams Kmodel_allply - predictors = c("discharge.daily"), + predictors = c("discharge.daily"), transforms = c(K600='log', date=NA, velocity.daily="log", discharge.daily="log"), other_args = c(), - - + + ## Sim - + # multi-day simulation parameters. already above for bayes: # K600_lnQ_nodes_centers, K600_lnQ_nodediffs_sdlog K600_lnQ_cnode_meanlog = log(6), # distrib for the y=K600 values of the middle (or just past middle) node @@ -472,7 +474,7 @@ specs <- function( sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) }, - + # daily simulation parameters discharge_daily = function(n, ...) rnorm(n, 20, 3), DO_mod_1 = NULL, @@ -482,29 +484,29 @@ specs <- function( alpha = function(n, ...) pmax(0, rnorm(n, 0.0001, 0.00002)), ER_daily = function(n, ...) pmin(0, rnorm(n, -10, 5)), ER20 = function(n, ...) pmin(0, rnorm(n, -10, 4)), - + # sub-daily simulation parameters err_obs_sigma = 0.01, err_obs_phi = 0, err_proc_sigma = 0.2, err_proc_phi = 0, err_round = NA, - + # simulation replicability sim_seed = NA - + ) { - + # make it easier to enter custom specs by creating the type-specific default if model_name %in% 'mle', etc. if(model_name %in% eval(formals(mm_name)$type)) model_name <- mm_name(type=model_name) - + # check the validity of the model_name against the list of officially accepted model names mm_validate_name(model_name) - + # parse the model_name features <- mm_parse_name(model_name, expand=TRUE) - + # collect info about the arguments required <- 'model_name' all_possible <- names(formals(specs)) @@ -512,11 +514,11 @@ specs <- function( yes_missing <- all_possible[!(all_possible %in% not_missing)] prefer_missing <- setdiff(all_possible[sapply(formals(specs), is.symbol)], 'params_out') # the arguments w/o defaults, mostly prefer_not_missing <- if(features$type == 'bayes' && features$GPP_fun == 'satlight') { - c('alpha_meanlog', 'alpha_sdlog', 'Pmax_mu', 'Pmax_sigma') + c('alpha_meanlog', 'alpha_sdlog', 'Pmax_mu', 'Pmax_sigma') } else { c() # could be made more extensive } - + # argument checks if(any(required %in% yes_missing)) stop("missing and required argument: ", paste(required[required %in% yes_missing], collapse=", ")) @@ -532,16 +534,16 @@ specs <- function( if(length(redundant) > 0) { warning("argument[s] that should usually be specified in revise() rather than specs(): ", paste(redundant, collapse=", ")) } - + # collect the defaults + directly specified arguments all_specs <- as.list(environment()) - + # copy/calculate arguments as appropriate to the model specs <- list() switch( features$type, 'bayes' = { - + # list the specs that will make it all the way to the Stan model as data all_specs$params_in <- c( switch( @@ -566,27 +568,27 @@ specs <- function( if(features$err_proc_iid) 'err_proc_iid_sigma_scale', if(features$err_proc_GPP) 'err_mult_GPP_sdlog_sigma' ) - + # list all needed arguments included <- c( # model setup 'model_name', 'engine', 'split_dates', 'keep_mcmcs', 'keep_mcmc_data', - + # date ply day_tests 'day_start', 'day_end', 'day_tests', 'required_timestep', - - # discharge binning parameters are not params_in, though they're + + # discharge binning parameters are not params_in, though they're # conceptually related and therefore colocated in formals(specs) if(features$pool_K600_type == 'binned') c('K600_lnQ_nodes_centers'), - + # params_in is both a vector of specs to include and a vector to include in specs all_specs$params_in, 'params_in', - + # inheritParams runstan_bayes - 'params_out', 'n_chains', 'n_cores', + 'params_out', 'n_chains', 'n_cores', 'burnin_steps', 'saved_steps', 'thin_steps', 'verbose' ) - + # compute some arguments if('engine' %in% yes_missing) { all_specs$engine <- features$engine @@ -595,7 +597,7 @@ specs <- function( all_specs$split_dates <- switch( features$pool_K600_type, 'none' = FALSE, # pretty sure FALSE is faster. also allows hierarchical error terms - 'normal'=, 'linear'=, 'binned' = FALSE, + 'normal'=, 'linear'=, 'binned' = FALSE, stop("unknown pool_K600; unsure how to set split_dates")) } if(features$pool_K600_type == 'binned') { @@ -616,7 +618,7 @@ specs <- function( none=c(), normal=c('K600_daily_predlog'), linear=c('K600_daily_predlog', 'lnK600_lnQ_intercept', 'lnK600_lnQ_slope'), - binned=c('K600_daily_predlog', 'lnK600_lnQ_nodes')), + binned=c('K600_daily_predlog', 'lnK600_lnQ_nodes')), if(features$pool_K600_sd == 'fitted') switch( features$pool_K600_type, @@ -627,10 +629,10 @@ specs <- function( if(features$err_proc_iid) c('err_proc_iid_sigma', 'err_proc_iid'), if(features$err_proc_GPP) c('err_proc_GPP', 'GPP_pseudo_R2')) } - + # check for errors/inconsistencies model_path <- tryCatch( - mm_locate_filename(model_name), + mm_locate_filename(model_name), error=function(e) { warning(e) return(model_name) @@ -639,15 +641,15 @@ specs <- function( stop('engine must be specified for Bayesian models') }, - 'bayes_2station' = { + 'bayes_2s' = { - # bayes_2station has a single fixed model structure (see + # bayes_2s has a single fixed model structure (see # inst/models/b2_np_oi_tr_plrckm.stan), so params_in/params_out are # hardcoded here rather than built up from pool_K600/GPP_fun/ER_fun/ # err_* toggles as in the 'bayes' case above # the six scalar prior hyperparameters spliced into the Stan data list - # by mm_ts_prep_data() + # by prepdata_bayes_2s() all_specs$params_in <- c( 'GPP_daily_mu', 'GPP_daily_sigma', 'ER_daily_mu', 'ER_daily_sigma', @@ -697,15 +699,15 @@ specs <- function( # determine which init values will be needed . <- '.dplyr.var' init.needs <- paste0('init.', get_param_names(model_name)$required) - + # list all needed arguments included <- c('model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep', init.needs) - }, + }, 'night' = { # list all needed arguments included <- c('model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep') - + # some different defaults for night relative to other models if('day_start' %in% yes_missing) { all_specs$day_start <- 12 @@ -716,18 +718,18 @@ specs <- function( if('day_tests' %in% yes_missing) { all_specs$day_tests <- c(day_tests, 'include_sunset') } - - }, + + }, 'Kmodel' = { # list all needed arguments included <- c( 'model_name', 'engine', 'day_start', 'day_end', 'day_tests', 'required_timestep', 'weights', 'filters', 'predictors', 'transforms', 'other_args') - + if('engine' %in% yes_missing) { all_specs$engine <- features$engine } - + # some different defaults for each engine, because no one set of defaults # makes sense for all engines #if('weights' %in% yes_missing) all_specs$weights <- c("K600/CI") # same for all, so use default as in Usage @@ -752,12 +754,12 @@ specs <- function( if('other_args' %in% yes_missing) all_specs$other_args <- list(possible_args=names(formals('loess'))[-which(names(formals('loess')) %in% c('formula','data','weights'))]) } ) - + }, 'sim' = { # determine which daily parameters will be needed par_needs <- gsub('\\.', '_', unlist(get_param_names(model_name)[c('optional','required')])) - + # list all needed arguments included <- c( 'model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep', @@ -766,23 +768,23 @@ specs <- function( none=c(), normal=stop("pool_K600='normal' unavailable for now; try 'binned' instead"), linear=stop("pool_K600='linear' unavailable for now; try 'binned' instead"), # 'discharge_daily', etc. - binned=c('K600_lnQ_nodes_centers', + binned=c('K600_lnQ_nodes_centers', 'K600_lnQ_cnode_meanlog', 'K600_lnQ_cnode_sdlog', 'K600_lnQ_nodediffs_meanlog', 'K600_lnQ_nodediffs_sdlog', 'lnK600_lnQ_nodes')), par_needs, 'err_round', 'sim_seed') - + if(features$pool_K600 == 'binned') { if('K600_lnQ_nodes_centers' %in% yes_missing) # override the default, which is for 'bayes' rather than 'sim' all_specs$K600_lnQ_nodes_centers <- function(discharge.daily, ...) calc_bins(log(discharge.daily), 'width', width=0.2)$bounds } } ) - + # stop if truly irrelevant arguments were given - if(length(irrelevant <- not_missing[!(not_missing %in% included)]) > 0) + if(length(irrelevant <- not_missing[!(not_missing %in% included)]) > 0) stop("irrelevant argument: ", paste(irrelevant, collapse=", ")) - + # return just the arguments we actually need add_specs_class(all_specs[included]) - + } diff --git a/man/specs.Rd b/man/specs.Rd index d0e64fa..58e6542 100644 --- a/man/specs.Rd +++ b/man/specs.Rd @@ -42,12 +42,11 @@ specs( K600_lnQ_nodediffs_sdlog = 0.5, K600_lnQ_nodes_meanlog = rep(log(12), length(K600_lnQ_nodes_centers)), K600_lnQ_nodes_sdlog = rep(1.32, length(K600_lnQ_nodes_centers)), - K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none = 1, - normal_sdfixed = 0.05, NA), + K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none = 1, normal_sdfixed + = 0.05, NA), K600_daily_sigma = switch(mm_parse_name(model_name)$pool_K600, linear_sdfixed = 10, binned_sdfixed = 5, NA), - K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal = 0.05, - NA), + K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal = 0.05, NA), K600_daily_sigma_sigma = switch(mm_parse_name(model_name)$pool_K600, linear = 1.2, binned = 0.24, NA), err_obs_iid_sigma_scale = 0.03, @@ -56,6 +55,8 @@ specs( err_proc_acor_phi_beta = 1, err_proc_acor_sigma_scale = 1, err_mult_GPP_sdlog_sigma = 1, + K600_lnorm_meanlog = 2.484907, + K600_lnorm_sdlog = 1, params_in, params_out, n_chains = 4, @@ -74,9 +75,11 @@ specs( K600_lnQ_cnode_sdlog = 1, K600_lnQ_nodediffs_meanlog = 0.2, lnK600_lnQ_nodes = function(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, - K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog, ...) { - sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, - K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) }, + K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog, ...) { + + sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, + K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) + }, discharge_daily = function(n, ...) rnorm(n, 20, 3), DO_mod_1 = NULL, K600_daily = function(n, K600_daily_predlog = log(10), ...) pmax(0, rnorm(n, @@ -337,6 +340,13 @@ estimate GPP_inst. The effect is a special kind of process error that is proportional to light (with noise) and is applied to GPP rather than to dDO/dt.} +\item{K600_lnorm_meanlog}{hyperparameter for \code{type='bayes_2s'}. +The mean of a lognormal prior distribution for K600_daily.} + +\item{K600_lnorm_sdlog}{hyperparameter for \code{type='bayes_2s'}. The +standard deviation parameter of a lognormal prior distribution for +K600_daily.} + \item{params_in}{Character vector of hyperparameters to pass from the specs list into the data list for the MCMC run. Will be automatically generated during the specs() call; need only be revised if you're using a custom From 5f3030e4616287977ea97f7c9b10bf12b883c499 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Tue, 28 Jul 2026 16:39:31 -0700 Subject: [PATCH 11/17] Extract shared matrix-pivot/contiguous-sort logic 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. --- DESCRIPTION | 1 + R/metab_bayes.R | 12 ++--- R/metab_bayes_2s.R | 16 +++---- R/mm_time_by_date_matrix.R | 43 +++++++++++++++++ man/metab_bayes_2s.Rd | 3 +- man/mm_check_dates_contiguous.Rd | 27 +++++++++++ man/mm_time_by_date_matrix.Rd | 28 ++++++++++++ tests/testthat/test-metab_bayes.R | 76 +++++++++++++++++++++++++++++++ 8 files changed, 188 insertions(+), 18 deletions(-) create mode 100644 R/mm_time_by_date_matrix.R create mode 100644 man/mm_check_dates_contiguous.Rd create mode 100644 man/mm_time_by_date_matrix.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 5607700..bb79d8d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -96,6 +96,7 @@ Collate: 'specs-class.R' 'metab_model-class.R' 'metab_Kmodel.R' + 'mm_time_by_date_matrix.R' 'metab_bayes.R' 'metab_bayes_2s.R' 'metab_inputs.R' diff --git a/R/metab_bayes.R b/R/metab_bayes.R index 9db1c0d..c941299 100644 --- a/R/metab_bayes.R +++ b/R/metab_bayes.R @@ -1,4 +1,4 @@ -#' @include metab_model-class.R +#' @include metab_model-class.R mm_time_by_date_matrix.R NULL #' Basic Bayesian metabolism model fitting function @@ -479,15 +479,13 @@ prepdata_bayes <- function( ) stop("dates have differing numbers of rows; observations cannot be combined in matrix") } - time_by_date_matrix <- function(vec) { - matrix(data=vec, nrow=num_daily_obs, ncol=num_dates, byrow=FALSE) - } + time_by_date_matrix <- mm_time_by_date_matrix(num_daily_obs, num_dates) # double-check that our dates are going to line up with the input dates. this # should be redundant w/ above date_table checks, so just being extra careful - obs_dates <- time_by_date_matrix(format(data$date, format="%Y-%m-%d")) - unique_dates <- apply(obs_dates, MARGIN=2, FUN=function(timevec) unique(timevec)) - if(!all.equal(unique_dates, names(date_table))) stop("couldn't fit given dates into matrix") + mm_check_dates_contiguous( + time_by_date_matrix(format(data$date, format="%Y-%m-%d")), date_table, + "couldn't fit given dates into matrix") # confirm that every day has the same modal timestep and put a value on that # timestep. the tolerance for uniqueness within each day is set by the default diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 94461ed..2537369 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -289,9 +289,9 @@ prepdata_bayes_2s <- function(data, specs=NULL) { travel_time = data$travel.time[keep] ) - # pivot into n_obs x n_days matrices, one column per unique date, following - # the same time_by_date_matrix approach used for the 1-station Stan models - # (see prepdata_bayes() in metab_bayes.R) + # pivot into n_obs x n_days matrices, one column per unique date, using the + # same mm_time_by_date_matrix()/mm_check_dates_contiguous() helpers shared + # with prepdata_bayes() (see mm_time_by_date_matrix.R) date_vec <- as.character(as.Date(modeled$solar.time)) date_table <- table(date_vec) n_days <- length(date_table) @@ -304,16 +304,14 @@ prepdata_bayes_2s <- function(data, specs=NULL) { } n_obs <- n_obs_per_day - to_matrix <- function(vec) matrix(vec, nrow=n_obs, ncol=n_days, byrow=FALSE) + to_matrix <- mm_time_by_date_matrix(n_obs, n_days) # confirm each date occupies a contiguous block of rows, i.e., that data # was sorted by solar.time; otherwise the matrix pivot below would silently # scramble which rows belong to which date - date_mat <- to_matrix(date_vec) - unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) - if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { - stop('data must be sorted by solar.time so that each date occupies a contiguous block of rows') - } + mm_check_dates_contiguous( + to_matrix(date_vec), date_table, + 'data must be sorted by solar.time so that each date occupies a contiguous block of rows') list( n_obs = n_obs, diff --git a/R/mm_time_by_date_matrix.R b/R/mm_time_by_date_matrix.R new file mode 100644 index 0000000..bbdbde3 --- /dev/null +++ b/R/mm_time_by_date_matrix.R @@ -0,0 +1,43 @@ +#' Build a closure that pivots a per-row vector into a time-by-date matrix +#' +#' Both \code{prepdata_bayes} (one-station) and \code{prepdata_bayes_2s} +#' (two-station) reshape several per-row vectors (DO, depth, light, etc.) +#' into an obs-per-day x num-days matrix, assuming that the input vector is +#' already sorted so that each date occupies a contiguous block of +#' \code{n_per_group} rows. This function returns the reshaping closure; +#' \code{\link{mm_check_dates_contiguous}} verifies the contiguous-block +#' assumption actually holds. +#' +#' @param n_per_group the number of rows per date (must be the same for +#' every date; callers are responsible for having already confirmed this) +#' @param n_groups the number of distinct dates +#' @return a function of one argument, \code{vec}, that reshapes \code{vec} +#' into an \code{n_per_group} x \code{n_groups} matrix, filling by column +#' @keywords internal +mm_time_by_date_matrix <- function(n_per_group, n_groups) { + function(vec) matrix(vec, nrow=n_per_group, ncol=n_groups, byrow=FALSE) +} + +#' Confirm that each date occupies a contiguous block of rows +#' +#' Used immediately after pivoting a date vector with the closure from +#' \code{\link{mm_time_by_date_matrix}}, to catch input that wasn't actually +#' sorted by date/time before pivoting (which would otherwise let the matrix +#' reshape silently scramble which rows belong to which date). Shared by +#' \code{prepdata_bayes} and \code{prepdata_bayes_2s}. +#' +#' @param date_mat the date-identifier vector (as used to build +#' \code{date_table}) already pivoted via +#' \code{\link{mm_time_by_date_matrix}}'s closure +#' @param date_table a table of date counts, as from \code{table(date_vec)}, +#' giving the expected date for each column of \code{date_mat} +#' @param error_message the message to pass to \code{stop()} if the dates +#' are not contiguous; left to the caller so each can keep its own wording +#' @keywords internal +mm_check_dates_contiguous <- function(date_mat, date_table, error_message) { + unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) + if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { + stop(error_message) + } + invisible(TRUE) +} diff --git a/man/metab_bayes_2s.Rd b/man/metab_bayes_2s.Rd index 4f7571e..bab4567 100644 --- a/man/metab_bayes_2s.Rd +++ b/man/metab_bayes_2s.Rd @@ -66,8 +66,7 @@ previous date's last rows. the previous day's light conditions from influencing the following day's metabolism estimate. There must also be enough lead-in observations of upstream DO before the first row of \code{data} to - cover the longest travel time in the dataset, given the (median) - timestep of \code{data$solar.time}. + cover the longest travel time in the dataset. } \seealso{ diff --git a/man/mm_check_dates_contiguous.Rd b/man/mm_check_dates_contiguous.Rd new file mode 100644 index 0000000..8f2821b --- /dev/null +++ b/man/mm_check_dates_contiguous.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_time_by_date_matrix.R +\name{mm_check_dates_contiguous} +\alias{mm_check_dates_contiguous} +\title{Confirm that each date occupies a contiguous block of rows} +\usage{ +mm_check_dates_contiguous(date_mat, date_table, error_message) +} +\arguments{ +\item{date_mat}{the date-identifier vector (as used to build +\code{date_table}) already pivoted via +\code{\link{mm_time_by_date_matrix}}'s closure} + +\item{date_table}{a table of date counts, as from \code{table(date_vec)}, +giving the expected date for each column of \code{date_mat}} + +\item{error_message}{the message to pass to \code{stop()} if the dates +are not contiguous; left to the caller so each can keep its own wording} +} +\description{ +Used immediately after pivoting a date vector with the closure from +\code{\link{mm_time_by_date_matrix}}, to catch input that wasn't actually +sorted by date/time before pivoting (which would otherwise let the matrix +reshape silently scramble which rows belong to which date). Shared by +\code{prepdata_bayes} and \code{prepdata_bayes_2s}. +} +\keyword{internal} diff --git a/man/mm_time_by_date_matrix.Rd b/man/mm_time_by_date_matrix.Rd new file mode 100644 index 0000000..762b5ce --- /dev/null +++ b/man/mm_time_by_date_matrix.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_time_by_date_matrix.R +\name{mm_time_by_date_matrix} +\alias{mm_time_by_date_matrix} +\title{Build a closure that pivots a per-row vector into a time-by-date matrix} +\usage{ +mm_time_by_date_matrix(n_per_group, n_groups) +} +\arguments{ +\item{n_per_group}{the number of rows per date (must be the same for +every date; callers are responsible for having already confirmed this)} + +\item{n_groups}{the number of distinct dates} +} +\value{ +a function of one argument, \code{vec}, that reshapes \code{vec} + into an \code{n_per_group} x \code{n_groups} matrix, filling by column +} +\description{ +Both \code{prepdata_bayes} (one-station) and \code{prepdata_bayes_2s} +(two-station) reshape several per-row vectors (DO, depth, light, etc.) +into an obs-per-day x num-days matrix, assuming that the input vector is +already sorted so that each date occupies a contiguous block of +\code{n_per_group} rows. This function returns the reshaping closure; +\code{\link{mm_check_dates_contiguous}} verifies the contiguous-block +assumption actually holds. +} +\keyword{internal} diff --git a/tests/testthat/test-metab_bayes.R b/tests/testthat/test-metab_bayes.R index 929872e..d5b527b 100644 --- a/tests/testthat/test-metab_bayes.R +++ b/tests/testthat/test-metab_bayes.R @@ -6,6 +6,82 @@ # skip_on_appveyor() # skip_if_not_installed('deSolve') + +# prepdata_bayes() --------------------------------------------------------- +# +# Direct, fast (no Stan) unit tests for prepdata_bayes()'s matrix-pivot and +# contiguous-sort-check behavior. Written to pin current behavior ahead of +# extracting shared logic with prepdata_bayes_2s() (see mm_time_by_date_matrix.R). + +# Build a minimal, traceable two-day data.frame: 4 rows/day, DO.obs == +# original row index so the pivoted matrix's values can be checked exactly. +make_bayes_prepdata <- function(n_per_day=4, n_days=2, shuffle=FALSE) { + n_total <- n_per_day * n_days + solar.time <- as.POSIXct("2050-06-01 00:00:00", tz="UTC") + + as.difftime( + rep((seq_len(n_per_day) - 1) * (24 / n_per_day), n_days) + + rep((seq_len(n_days) - 1) * 24, each=n_per_day), + units="hours") + dat <- data.frame( + solar.time = solar.time, + date = as.Date(solar.time), + DO.obs = seq_len(n_total), + DO.sat = rep(10, n_total), + depth = rep(1, n_total), + temp.water = rep(20, n_total), + light = rep(100, n_total) + ) + if(shuffle) { + # interleave day1/day2 rows to break contiguity while keeping an equal + # row count per date (so the earlier date_table-based check still passes + # and only the contiguous-sort check is exercised) + dat <- dat[order(rep(seq_len(n_per_day), n_days)), ] + } + dat +} + +test_that("prepdata_bayes() pivots data into the expected n x d matrix (dims and values)", { + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2) + + out <- prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp) + + expect_equal(out$d, 2) + expect_equal(out$n, 4) + expect_equal(dim(out$DO_obs), c(4, 2)) + expect_equal(dim(out$DO_sat), c(4, 2)) + expect_equal(dim(out$depth), c(4, 2)) + # values: matrix(vec, nrow=4, ncol=2, byrow=FALSE) fills column-wise, so + # day 1 = original rows 1:4, day 2 = original rows 5:8 + expect_equal(out$DO_obs[,1], as.numeric(1:4)) + expect_equal(out$DO_obs[,2], as.numeric(5:8)) +}) + +test_that("prepdata_bayes()'s contiguous-sort check passes for properly sorted data", { + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2) + + expect_silent(prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp)) +}) + +test_that("prepdata_bayes() errors for data that isn't sorted by date", { + # HISTORY: prior to the mm_time_by_date_matrix()/mm_check_dates_contiguous() + # extraction, this check was `if(!all.equal(unique_dates, names(date_table))) + # stop("couldn't fit given dates into matrix")`, unguarded by isTRUE(). + # Confirmed empirically pre-refactor: whenever the dates were truly + # non-contiguous, all.equal() returned a character vector (not TRUE), and + # `!` on a character vector always errors with "invalid argument type" in + # base R -- so the intended "couldn't fit given dates into matrix" message + # was actually unreachable. The shared mm_check_dates_contiguous() helper + # fixes this as a side effect of the extraction, so this test now + # asserts the originally-intended message. + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2, shuffle=TRUE) + + expect_error(prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp), "couldn't fit given dates into matrix") +}) + + manual_test4 <- function() { library(streamMetabolizer) From 4103b713ee2fb6a3627ca176ab921a34d3e9365f Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Thu, 30 Jul 2026 14:50:02 -0700 Subject: [PATCH 12/17] Extract shared core-count logic into mm_determine_cores() 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() --- DESCRIPTION | 1 + R/metab_bayes.R | 5 +-- R/metab_bayes_2s.R | 4 +-- R/mm_determine_cores.R | 27 ++++++++++++++ man/mm_determine_cores.Rd | 30 ++++++++++++++++ tests/testthat/test-mm_determine_cores.R | 46 ++++++++++++++++++++++++ 6 files changed, 106 insertions(+), 7 deletions(-) create mode 100644 R/mm_determine_cores.R create mode 100644 man/mm_determine_cores.Rd create mode 100644 tests/testthat/test-mm_determine_cores.R diff --git a/DESCRIPTION b/DESCRIPTION index bb79d8d..b3c27f8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -110,6 +110,7 @@ Collate: 'metab_sim.R' 'mm_check_mcmc_file.R' 'mm_data.R' + 'mm_determine_cores.R' 'mm_filter_dates.R' 'mm_filter_hours.R' 'mm_filter_valid_days.R' diff --git a/R/metab_bayes.R b/R/metab_bayes.R index c941299..d0dcac0 100644 --- a/R/metab_bayes.R +++ b/R/metab_bayes.R @@ -623,10 +623,7 @@ runstan_bayes <- function( verbose=FALSE, ...) { # determine how many cores to use - tot_cores <- detectCores() - if (!is.finite(tot_cores)) { tot_cores <- 1 } - n_cores <- min(tot_cores, n_cores) - if(verbose) message(paste0("MCMC (","Stan","): requesting ",n_chains," chains on ",n_cores," of ",tot_cores," available cores")) + n_cores <- mm_determine_cores(n_cores, n_chains=n_chains, verbose=verbose) # stan() can't find its own function cpp_object_initializer() unless the # namespace is loaded. requireNamespace is somehow not doing this. Thoughts diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 2537369..0000983 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -101,9 +101,7 @@ metab_bayes_2s <- function( specs$model_path <- mm_locate_filename(specs$model_name) # determine how many cores to use, as in runstan_bayes() - tot_cores <- parallel::detectCores() - if(!is.finite(tot_cores)) tot_cores <- 1 - n_cores <- min(tot_cores, specs$n_cores) + n_cores <- mm_determine_cores(specs$n_cores, n_chains=specs$n_chains, verbose=specs$verbose) # Fit the model, collecting errors/warnings as strings rather than # letting a bad dataset halt execution without reporting anything back diff --git a/R/mm_determine_cores.R b/R/mm_determine_cores.R new file mode 100644 index 0000000..7c9ea33 --- /dev/null +++ b/R/mm_determine_cores.R @@ -0,0 +1,27 @@ +#' Determine how many cores to use for an MCMC run +#' +#' Shared by \code{runstan_bayes} (one-station) and \code{metab_bayes_2s} +#' (two-station): detects the number of cores available on the machine, +#' falls back to 1 if detection fails, and caps the requested core count at +#' whatever is actually available. +#' +#' @param n_cores the number of cores requested for this run +#' @param n_chains the number of chains being requested, used only to +#' reconstruct the verbose status message; if NULL (the default), the +#' chains clause is omitted from the message. Ignored when +#' \code{verbose=FALSE} +#' @param verbose logical. if TRUE, emit a status message reporting the +#' number of cores requested vs. available +#' @return the number of cores to actually use, i.e. +#' \code{min(detected_cores, n_cores)} +#' @keywords internal +mm_determine_cores <- function(n_cores, n_chains=NULL, verbose=FALSE) { + tot_cores <- parallel::detectCores() + if(!is.finite(tot_cores)) tot_cores <- 1 + n_cores <- min(tot_cores, n_cores) + if(verbose) { + chains_clause <- if(is.null(n_chains)) "" else paste0(n_chains," chains on ") + message(paste0("MCMC (","Stan","): requesting ",chains_clause,n_cores," of ",tot_cores," available cores")) + } + n_cores +} diff --git a/man/mm_determine_cores.Rd b/man/mm_determine_cores.Rd new file mode 100644 index 0000000..cedda66 --- /dev/null +++ b/man/mm_determine_cores.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_determine_cores.R +\name{mm_determine_cores} +\alias{mm_determine_cores} +\title{Determine how many cores to use for an MCMC run} +\usage{ +mm_determine_cores(n_cores, n_chains = NULL, verbose = FALSE) +} +\arguments{ +\item{n_cores}{the number of cores requested for this run} + +\item{n_chains}{the number of chains being requested, used only to +reconstruct the verbose status message; if NULL (the default), the +chains clause is omitted from the message. Ignored when +\code{verbose=FALSE}} + +\item{verbose}{logical. if TRUE, emit a status message reporting the +number of cores requested vs. available} +} +\value{ +the number of cores to actually use, i.e. + \code{min(detected_cores, n_cores)} +} +\description{ +Shared by \code{runstan_bayes} (one-station) and \code{metab_bayes_2s} +(two-station): detects the number of cores available on the machine, +falls back to 1 if detection fails, and caps the requested core count at +whatever is actually available. +} +\keyword{internal} diff --git a/tests/testthat/test-mm_determine_cores.R b/tests/testthat/test-mm_determine_cores.R new file mode 100644 index 0000000..4c8c615 --- /dev/null +++ b/tests/testthat/test-mm_determine_cores.R @@ -0,0 +1,46 @@ +context("mm_determine_cores") + +test_that("falls back to 1 core when detectCores() is non-finite", { + local_mocked_bindings(detectCores = function() NA_real_, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=5), 1) + + local_mocked_bindings(detectCores = function() Inf, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=5), 1) +}) + +test_that("caps the requested core count at the number detected", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=10), 4) + expect_equal(mm_determine_cores(n_cores=2), 2) + expect_equal(mm_determine_cores(n_cores=4), 4) +}) + +test_that("emits a status message only when verbose=TRUE", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + + expect_message( + result <- mm_determine_cores(n_cores=10, n_chains=3, verbose=TRUE), + "requesting 3 chains on 4 of 4 available cores", + fixed=TRUE) + expect_equal(result, 4) + + expect_no_message(mm_determine_cores(n_cores=10, n_chains=3, verbose=FALSE)) +}) + +test_that("verbose message matches runstan_bayes()'s original wording exactly", { + local_mocked_bindings(detectCores = function() 8, .package = "parallel") + + expect_message( + mm_determine_cores(n_cores=4, n_chains=4, verbose=TRUE), + "MCMC (Stan): requesting 4 chains on 4 of 8 available cores", + fixed=TRUE) +}) + +test_that("omits the chains clause when n_chains is not supplied", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + + expect_message( + mm_determine_cores(n_cores=10, verbose=TRUE), + "MCMC (Stan): requesting 4 of 4 available cores", + fixed=TRUE) +}) From 13b8c8f94f1fdaf9f3e6d758666f6ac3b65f27eb Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Sat, 1 Aug 2026 08:53:27 -0700 Subject: [PATCH 13/17] Align failed-Stan-run behavior in metab_bayes_2s() with one-station 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). --- R/metab_bayes_2s.R | 61 +++++++++++++++------------- tests/testthat/test-metab_bayes_2s.R | 32 +++++++++++++++ 2 files changed, 64 insertions(+), 29 deletions(-) diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 0000983..2828039 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -122,37 +122,40 @@ metab_bayes_2s <- function( split=specs$verbose) if(stanfit@mode == 2L) { - stop(paste(utils::capture.output(print(stanfit)), collapse='\n')) + # mirror runstan_bayes()'s warn-and-continue pattern for a failed + # run: report the diagnostic as a warning (caught below) and skip + # the post-processing steps, which all assume a successful fit + warning(paste(utils::capture.output(print(stanfit)), collapse='\n')) + } else { + # format the Stan summary matrix into per-variable data.frames + stan_mat <- rstan::summary(stanfit)$summary + mcmc_out <- format_mcmc_mat_nosplit( + stan_mat, data_list$n_days, data_list$n_obs, specs$model_name, + keep_mcmc=isTRUE(specs$keep_mcmcs), stanfit) + + # daily GPP/ER/K600 estimates: join Stan's date_index back to dates + date_index <- time_index <- index <- '.dplyr.var' + daily <- mcmc_out$daily %>% + dplyr::left_join(date_df, by='date_index') %>% + dplyr::select(-date_index, -time_index, -index) %>% + dplyr::select(date, dplyr::everything()) + + # instantaneous DO.mod.down estimates come from the 'metab' Stan + # transformed parameter (posterior median), which format_mcmc_mat_nosplit() + # buckets by row count rather than by name since 'metab' isn't in its + # par_homes lookup table; find that bucket by its column names instead + is_metab_bucket <- sapply(mcmc_out, function(df) is.data.frame(df) && any(grepl('^metab_', names(df)))) + metab_bucket_name <- names(mcmc_out)[is_metab_bucket][1] + if(is.na(metab_bucket_name)) { + stop("could not find 'metab' in the Stan output; check that specs$params_out includes 'metab'") + } + inst <- mcmc_out[[metab_bucket_name]] %>% + dplyr::select(date_index, time_index, DO.mod.down=metab_50pct) %>% + dplyr::inner_join(obs_index_df, by=c('date_index','time_index')) %>% + dplyr::select(solar.time, DO.obs.down, DO.mod.down) %>% + dplyr::arrange(solar.time) } - # format the Stan summary matrix into per-variable data.frames - stan_mat <- rstan::summary(stanfit)$summary - mcmc_out <- format_mcmc_mat_nosplit( - stan_mat, data_list$n_days, data_list$n_obs, specs$model_name, - keep_mcmc=isTRUE(specs$keep_mcmcs), stanfit) - - # daily GPP/ER/K600 estimates: join Stan's date_index back to dates - date_index <- time_index <- index <- '.dplyr.var' - daily <- mcmc_out$daily %>% - dplyr::left_join(date_df, by='date_index') %>% - dplyr::select(-date_index, -time_index, -index) %>% - dplyr::select(date, dplyr::everything()) - - # instantaneous DO.mod.down estimates come from the 'metab' Stan - # transformed parameter (posterior median), which format_mcmc_mat_nosplit() - # buckets by row count rather than by name since 'metab' isn't in its - # par_homes lookup table; find that bucket by its column names instead - is_metab_bucket <- sapply(mcmc_out, function(df) is.data.frame(df) && any(grepl('^metab_', names(df)))) - metab_bucket_name <- names(mcmc_out)[is_metab_bucket][1] - if(is.na(metab_bucket_name)) { - stop("could not find 'metab' in the Stan output; check that specs$params_out includes 'metab'") - } - inst <- mcmc_out[[metab_bucket_name]] %>% - dplyr::select(date_index, time_index, DO.mod.down=metab_50pct) %>% - dplyr::inner_join(obs_index_df, by=c('date_index','time_index')) %>% - dplyr::select(solar.time, DO.obs.down, DO.mod.down) %>% - dplyr::arrange(solar.time) - }, error=function(err) { stop_strs <<- c(stop_strs, err$message) }), warning=function(war) { diff --git a/tests/testthat/test-metab_bayes_2s.R b/tests/testthat/test-metab_bayes_2s.R index 1c8be70..34abe0f 100644 --- a/tests/testthat/test-metab_bayes_2s.R +++ b/tests/testthat/test-metab_bayes_2s.R @@ -241,3 +241,35 @@ test_that("metab() fits a two-station model and predict_metab()/predict_DO() wor expect_s3_class(pdo, 'data.frame') expect_true(all(c('DO.obs.down','DO.mod.down') %in% names(pdo))) }) + +test_that("a failed Stan run (mode==2L) warns and continues, matching runstan_bayes()'s pattern, rather than erroring out", { + skip_if_not_installed('rstan') + + # stand in for rstan::stan()'s return value on a failed run: only the + # 'mode' slot is inspected by metab_bayes_2s() before deciding to skip + # post-processing, so a minimal S4 object with that slot is sufficient + setClass('fake_failed_stanfit', representation(mode='integer')) + fake_stanfit <- methods::new('fake_failed_stanfit', mode=2L) + testthat::local_mocked_bindings(stan=function(...) fake_stanfit, .package='rstan') + + dat <- make_ts_data() + sp <- specs( + mm_name('bayes_2s'), + n_chains=1, n_cores=1, burnin_steps=10, saved_steps=10, verbose=FALSE) + + expect_warning( + mm <- metab_bayes_2s(specs=sp, data=dat), + 'Modeling failed') + + expect_s4_class(mm, 'metab_bayes_2s') + fit <- mm@fit + expect_true(nrow(fit$daily) > 0) + expect_true(all(is.na(fit$daily$GPP_daily_50pct))) + expect_true(all(is.na(fit$daily$ER_daily_50pct))) + expect_true(all(is.na(fit$daily$K600_daily_50pct))) + expect_true(all(fit$daily$valid_day)) + expect_null(fit$inst) + expect_equal(length(fit$errors), 0) + expect_true(length(fit$warnings) > 0) + expect_true(any(grepl('fake_failed_stanfit', fit$warnings))) +}) From d270bfef070dcf9790a770e5fb41f6a8734a2aea Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Sun, 2 Aug 2026 09:43:17 -0700 Subject: [PATCH 14/17] Fix units notation to unitted-style 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. --- R/data.R | 8 ++++---- man/two_station_example.Rd | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/R/data.R b/R/data.R index b0a00ab..64c1c25 100644 --- a/R/data.R +++ b/R/data.R @@ -16,13 +16,13 @@ #' \describe{ #' \item{solar.time}{POSIXct timestamp, UTC} #' \item{DO.obs.up}{dissolved oxygen observed at the upstream station, -#' mg O2 / L} +#' mgO2 L^-1} #' \item{DO.sat.up}{dissolved oxygen at equilibrium saturation at the -#' upstream station, mg O2 / L} +#' upstream station, mgO2 L^-1} #' \item{DO.obs.down}{dissolved oxygen observed at the downstream -#' station, mg O2 / L} +#' station, mgO2 L^-1} #' \item{DO.sat.down}{dissolved oxygen at equilibrium saturation at the -#' downstream station, mg O2 / L} +#' downstream station, mgO2 L^-1} #' \item{light}{photosynthetically active radiation, umol m^-2 s^-1} #' \item{depth}{reach depth, m} #' \item{temp.water}{water temperature at the downstream station, degC} diff --git a/man/two_station_example.Rd b/man/two_station_example.Rd index da4ab53..7bcafb2 100644 --- a/man/two_station_example.Rd +++ b/man/two_station_example.Rd @@ -11,13 +11,13 @@ A data.frame with 2904 rows and the 9 columns expected by \describe{ \item{solar.time}{POSIXct timestamp, UTC} \item{DO.obs.up}{dissolved oxygen observed at the upstream station, - mg O2 / L} + mgO2 L^-1} \item{DO.sat.up}{dissolved oxygen at equilibrium saturation at the - upstream station, mg O2 / L} + upstream station, mgO2 L^-1} \item{DO.obs.down}{dissolved oxygen observed at the downstream - station, mg O2 / L} + station, mgO2 L^-1} \item{DO.sat.down}{dissolved oxygen at equilibrium saturation at the - downstream station, mg O2 / L} + downstream station, mgO2 L^-1} \item{light}{photosynthetically active radiation, umol m^-2 s^-1} \item{depth}{reach depth, m} \item{temp.water}{water temperature at the downstream station, degC} From 8a81e49487ab5b7c5f8e3da8ca6a5de26e0ba8b4 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Mon, 3 Aug 2026 08:05:31 -0700 Subject: [PATCH 15/17] Fix stale bayes_allply() data_all doc @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. --- R/metab_bayes.R | 8 +++++--- man/bayes_allply.Rd | 8 +++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/R/metab_bayes.R b/R/metab_bayes.R index d0dcac0..43d49c3 100644 --- a/R/metab_bayes.R +++ b/R/metab_bayes.R @@ -330,9 +330,11 @@ bayes_1ply <- function( #' Called from metab_bayes(). #' #' @param data_all data.frame of the form \code{mm_data(solar.time, DO.obs, -#' DO.sat, depth, temp.water, light)} and containing data for just one -#' estimation-day (this may be >24 hours but only yields estimates for one -#' 24-hour period) +#' DO.sat, depth, temp.water, light)} containing the full (possibly +#' multi-day) filtered dataset for this model - unlike \code{bayes_1ply()}'s +#' \code{data_ply}, which receives one estimation-day at a time, +#' \code{bayes_allply()} is called once with all valid dates together (used +#' when \code{specs$split_dates==FALSE}) #' @param data_daily_all data.frame of daily priors, if appropriate to the given #' model_path #' @param removed data.frame of dates that were removed and why diff --git a/man/bayes_allply.Rd b/man/bayes_allply.Rd index 46e5c25..c1c77e8 100644 --- a/man/bayes_allply.Rd +++ b/man/bayes_allply.Rd @@ -9,9 +9,11 @@ bayes_allply(data_all, data_daily_all, removed, specs) } \arguments{ \item{data_all}{data.frame of the form \code{mm_data(solar.time, DO.obs, -DO.sat, depth, temp.water, light)} and containing data for just one -estimation-day (this may be >24 hours but only yields estimates for one -24-hour period)} +DO.sat, depth, temp.water, light)} containing the full (possibly +multi-day) filtered dataset for this model - unlike \code{bayes_1ply()}'s +\code{data_ply}, which receives one estimation-day at a time, +\code{bayes_allply()} is called once with all valid dates together (used +when \code{specs$split_dates==FALSE})} \item{data_daily_all}{data.frame of daily priors, if appropriate to the given model_path} From e56233882d521271d5ceff410c6f89dd41289021 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Fri, 28 Aug 2026 09:30:13 -0700 Subject: [PATCH 16/17] Add attach.units deprecation support to predict_metab.metab_bayes_2s() 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. --- R/metab_bayes_2s.R | 17 ++++++++++++++++- R/mm_data.R | 3 +++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 2828039..1af01a9 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -415,9 +415,19 @@ get_params.metab_bayes_2s <- function( #' the two-station Stan model results. #' @export #' @import dplyr -predict_metab.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ...) { +#' @importFrom lifecycle deprecated is_present +#' @importFrom unitted get_units u +predict_metab.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ..., attach.units=deprecated()) { Var1 <- Var2 <- '.dplyr.var' + + # check units-related arguments + if (lifecycle::is_present(attach.units)) { + unitted_deprecate_warn("predict_metab(attach.units)") + } else { + attach.units <- FALSE + } + fit.names <- expand.grid(c('50pct','2.5pct','97.5pct'), c('GPP_daily','ER_daily','K600_daily'), stringsAsFactors=FALSE) %>% select(Var2, Var1) %>% apply(MARGIN=1, FUN=function(row) do.call(paste, c(as.list(row), list(sep='_')))) @@ -449,6 +459,11 @@ predict_metab.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA warnings=if(length(metab_model@fit$errors) > 0) NA else '', errors=if(length(metab_model@fit$errors) > 0) NA else '') + # attach.units if requested + if(attach.units) { + pred.units <- get_units(mm_data())[sapply(names(preds), function(x) strsplit(x, '\\.')[[1]][1], USE.NAMES=FALSE)] + preds <- u(preds, pred.units) + } preds } diff --git a/R/mm_data.R b/R/mm_data.R index e6009c4..da229d6 100644 --- a/R/mm_data.R +++ b/R/mm_data.R @@ -153,6 +153,9 @@ mm_data <- function(..., optional='none') { ER = u(-5,"gO2 m^-2 d^-1"), ER.lower = u(-6,"gO2 m^-2 d^-1"), ER.upper = u(-4,"gO2 m^-2 d^-1"), + K600 = u(10,"d^-1"), + K600.lower = u(4.5,"d^-1"), + K600.upper = u(15.6,"d^-1"), D = u(5,"gO2 m^-3 d^-1"), D.lower = u(5,"gO2 m^-3 d^-1"), D.upper = u(5,"gO2 m^-3 d^-1") From 4326791e310e09020fe1704f7b2b91d2384b3d80 Mon Sep 17 00:00:00 2001 From: Mike Dodrill Date: Fri, 28 Aug 2026 10:27:05 -0700 Subject: [PATCH 17/17] Doc-trim pass on branch's own roxygen (fix/vfts-stan-syntax) 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. --- R/metab_bayes_2s.R | 16 ++++++++-------- R/mm_validate_data.R | 11 +++++------ man/get_params.Rd | 12 ++++-------- man/mm_data.Rd | 18 ++++++++++++++++++ man/mm_parse_name.Rd | 10 +++++----- man/mm_validate_data_2station.Rd | 11 +++++------ man/predict_metab.Rd | 8 +++++++- man/streamMetabolizer.Rd | 20 ++++++++++++++++++++ 8 files changed, 72 insertions(+), 34 deletions(-) diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R index 1af01a9..d158270 100644 --- a/R/metab_bayes_2s.R +++ b/R/metab_bayes_2s.R @@ -351,14 +351,10 @@ setClass("metab_bayes_2s", contains="metab_bayes") #' @describeIn get_params Does the same Stan-output-to-streamMetabolizer -#' renaming as \code{get_params.metab_bayes}, but (unlike that method) does -#' not delegate the rest of the work to \code{get_params.metab_model} via -#' \code{NextMethod()}: that generic implementation looks up parameter -#' names via \code{get_param_names()}, which builds an ODE-based dDOdt -#' function using streamMetabolizer's one-station instantaneous-rate -#' framework (\code{ode_method}/\code{GPP_fun}/\code{ER_fun}/ -#' \code{deficit_src}) -- machinery that doesn't apply to the two-station -#' steady-state model's daily GPP/ER/K600 parameters. \code{fixed} +#' renaming as \code{get_params.metab_bayes}, but does not delegate to +#' \code{get_params.metab_model} via \code{NextMethod()}: the two-station +#' steady-state model's daily GPP/ER/K600 parameters don't fit that +#' generic's one-station, ODE-based parameter-name lookup. \code{fixed} #' column/star annotations (relevant only to models that can take fixed #' daily parameters from \code{data_daily}) are not supported here. #' @export @@ -366,6 +362,10 @@ setClass("metab_bayes_2s", contains="metab_bayes") get_params.metab_bayes_2s <- function( metab_model, date_start=NA, date_end=NA, uncertainty=c('sd','ci','none'), messages=TRUE, ...) { + # not delegated to get_params.metab_model via NextMethod(): that generic's + # parameter-name lookup builds an ODE-based dDOdt function from one-station's + # ode_method/GPP_fun/ER_fun/deficit_src specs, which don't exist for this + # steady-state model uncertainty <- match.arg(uncertainty) fit <- metab_model@fit$daily diff --git a/R/mm_validate_data.R b/R/mm_validate_data.R index 0b7bf90..5e6f14f 100644 --- a/R/mm_validate_data.R +++ b/R/mm_validate_data.R @@ -123,12 +123,11 @@ mm_validate_data <- function( #' Two-station-specific data validation #' -#' Checks the lead-in coverage requirement specific to -#' \code{\link{metab_bayes_2s}}: there must be enough lead-in observations of -#' upstream DO before the first modeled row to cover the longest travel time -#' in the dataset, given the (median) timestep of \code{data$solar.time}. -#' Column presence, timestamp validity, and travel.time bounds are expected -#' to have already been checked by \code{\link{mm_validate_data}}. +#' Checks the lead-in coverage requirement described in +#' \code{\link{metab_bayes_2s}}, using the (median) timestep of +#' \code{data$solar.time} to compute the required lag. Column presence, +#' timestamp validity, and travel.time bounds are expected to have already +#' been checked by \code{\link{mm_validate_data}}. #' #' @param data data.frame as returned by \code{\link{mm_validate_data}} for #' \code{\link{metab_bayes_2s}}: must contain \code{solar.time} and diff --git a/man/get_params.Rd b/man/get_params.Rd index 2a5b3c8..69526bd 100644 --- a/man/get_params.Rd +++ b/man/get_params.Rd @@ -129,14 +129,10 @@ to streamMetabolizer parameter names; otherwise the same as \code{get_params.metab_model} \item \code{get_params(metab_bayes_2s)}: Does the same Stan-output-to-streamMetabolizer -renaming as \code{get_params.metab_bayes}, but (unlike that method) does -not delegate the rest of the work to \code{get_params.metab_model} via -\code{NextMethod()}: that generic implementation looks up parameter -names via \code{get_param_names()}, which builds an ODE-based dDOdt -function using streamMetabolizer's one-station instantaneous-rate -framework (\code{ode_method}/\code{GPP_fun}/\code{ER_fun}/ -\code{deficit_src}) -- machinery that doesn't apply to the two-station -steady-state model's daily GPP/ER/K600 parameters. \code{fixed} +renaming as \code{get_params.metab_bayes}, but does not delegate to +\code{get_params.metab_model} via \code{NextMethod()}: the two-station +steady-state model's daily GPP/ER/K600 parameters don't fit that +generic's one-station, ODE-based parameter-name lookup. \code{fixed} column/star annotations (relevant only to models that can take fixed daily parameters from \code{data_daily}) are not supported here. diff --git a/man/mm_data.Rd b/man/mm_data.Rd index 7fe0459..af079a9 100644 --- a/man/mm_data.Rd +++ b/man/mm_data.Rd @@ -43,6 +43,24 @@ Produces a unitted data.frame with the column names, units, and equilibrium saturation \eqn{mg O[2] L^{-1}}{mg O2 / L}. Calculate using \link{calc_DO_sat}} + \item{ \code{DO.obs.up} dissolved oxygen concentration observations at the + upstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 / L}} + + \item{ \code{DO.sat.up} dissolved oxygen concentrations at equilibrium + saturation at the upstream station of a two-station reach, \eqn{mg O[2] + L^{-1}}{mg O2 / L}} + + \item{ \code{DO.obs.down} dissolved oxygen concentration observations at + the downstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 + / L}} + + \item{ \code{DO.sat.down} dissolved oxygen concentrations at equilibrium + saturation at the downstream station of a two-station reach, \eqn{mg O[2] + L^{-1}}{mg O2 / L}} + + \item{ \code{travel.time} reach travel time between the upstream and + downstream stations of a two-station reach, in days, \eqn{d}{d}} + \item{ \code{depth} stream depth, \eqn{m}{m}}. \item{ \code{temp.water} water temperature, \eqn{degC}}. diff --git a/man/mm_parse_name.Rd b/man/mm_parse_name.Rd index 670b22f..cc7a3f5 100644 --- a/man/mm_parse_name.Rd +++ b/man/mm_parse_name.Rd @@ -11,17 +11,17 @@ mm_parse_name(model_name, expand = FALSE) \item{expand}{logical: should additional columns such as model_name and pool_K600_type be added? If expand=TRUE then the result cannot be passed -directly back into mm_name, but the additional columns may be helpful for +directly back into mm_name, but the additional columns may be helpful for interpreting the model structure.} } \description{ -Returns a data.frame with one column per model structure detail and one row -per `model_name` supplied to this function. See \code{?\link{mm_name}} for a +Returns a data.frame with one column per model structure detail and one row +per `model_name` supplied to this function. See \code{?\link{mm_name}} for a description of each of the data.frame columns that is returned. } \details{ -Custom model files (for MCMC) may have additional characters after an -underscore at the end of the name and before the prefix. For example, +Custom model files (for MCMC) may have additional characters after an +underscore at the end of the name and before the prefix. For example, 'b_np_pcpi_eu_ko.stan' and 'b_np_pcpi_eu_ko_v2.stan' are parsed the same; the _v2 is ignored by this function. } diff --git a/man/mm_validate_data_2station.Rd b/man/mm_validate_data_2station.Rd index 61984e3..9b6c037 100644 --- a/man/mm_validate_data_2station.Rd +++ b/man/mm_validate_data_2station.Rd @@ -12,11 +12,10 @@ mm_validate_data_2station(data) \code{travel.time}, sorted ascending by \code{solar.time}.} } \description{ -Checks the lead-in coverage requirement specific to -\code{\link{metab_bayes_2s}}: there must be enough lead-in observations of -upstream DO before the first modeled row to cover the longest travel time -in the dataset, given the (median) timestep of \code{data$solar.time}. -Column presence, timestamp validity, and travel.time bounds are expected -to have already been checked by \code{\link{mm_validate_data}}. +Checks the lead-in coverage requirement described in +\code{\link{metab_bayes_2s}}, using the (median) timestep of +\code{data$solar.time} to compute the required lag. Column presence, +timestamp validity, and travel.time bounds are expected to have already +been checked by \code{\link{mm_validate_data}}. } \keyword{internal} diff --git a/man/predict_metab.Rd b/man/predict_metab.Rd index f53bf06..e4a5df9 100644 --- a/man/predict_metab.Rd +++ b/man/predict_metab.Rd @@ -27,7 +27,13 @@ predict_metab( attach.units = deprecated() ) -\method{predict_metab}{metab_bayes_2s}(metab_model, date_start = NA, date_end = NA, ...) +\method{predict_metab}{metab_bayes_2s}( + metab_model, + date_start = NA, + date_end = NA, + ..., + attach.units = deprecated() +) \method{predict_metab}{metab_model}( metab_model, diff --git a/man/streamMetabolizer.Rd b/man/streamMetabolizer.Rd index 25a0bf0..dd33051 100644 --- a/man/streamMetabolizer.Rd +++ b/man/streamMetabolizer.Rd @@ -106,3 +106,23 @@ See http://usgs-r.github.io/streamMetabolizer for vignettes on the web. } } +\seealso{ +Useful links: +\itemize{ + \item \url{https://github.com/USGS-R/streamMetabolizer} + \item \url{http://usgs-r.github.io/streamMetabolizer/} + \item Report bugs at \url{https://github.com/USGS-R/streamMetabolizer/issues} +} + +} +\author{ +\strong{Maintainer}: Alison P. Appling \email{aappling@usgs.gov} + +Authors: +\itemize{ + \item Robert O. Hall + \item Maite Arroita + \item Charles B. Yackulic +} + +}