Clean up workflow ergonomics, defer flat_prior-only inits, drop sero for now

.Rbuildignore

@@ -6,3 +6,4 @@
 ^README\.Rmd$
 ^\.positai$
 ^\.claude$
+^dev$

.gitignore

@@ -9,3 +9,4 @@ inst/doc
 figures/
 simdata/
 .positai
+README.html

NAMESPACE

@@ -1,26 +1,21 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method(as_matrix,epiwave_greta_timeseries)
 S3method(as_matrix,epiwave_timeseries)
+S3method(as_matrix,greta_timeseries)
 S3method(as_matrix,numeric)
+export(as_greta_timeseries)
 export(as_matrix)
 export(compute_reff)
 export(create_dow_priors)
-export(create_epiwave_fixed_timeseries)
-export(create_epiwave_greta_timeseries)
-export(create_epiwave_timeseries)
 export(create_infection_timeseries)
 export(create_observation_model)
-export(create_small_sero_model)
 export(define_observation_data)
 export(define_observation_model)
 export(estimate_reff)
 export(fit_waves)
 export(implement_day_of_week)
-export(inits_by_jurisdiction)
 export(new_convolution_matrix)
 export(plot_infection_traj)
-export(prepare_observation_data)
 export(stack_jurisdictions)
 importFrom(cowplot,panel_border)
 importFrom(cowplot,theme_cowplot)
@@ -35,7 +30,6 @@ importFrom(ggplot2,theme)
 importFrom(greta,"%*%")
 importFrom(greta,apply)
 importFrom(greta,as_data)
-importFrom(greta,binomial)
 importFrom(greta,dirichlet)
 importFrom(greta,lognormal)
 importFrom(greta,negative_binomial)

R/as_matrix.R

@@ -51,11 +51,11 @@ as_matrix.epiwave_timeseries <- function (data, target_infection_dates, ...) {
 }
 
 #' @export
-as_matrix.epiwave_greta_timeseries <- function (data, target_infection_dates, ...) {
+as_matrix.greta_timeseries <- function (data, target_infection_dates, ...) {
 
   if (!identical(as.Date(data$timeseries$date), as.Date(target_infection_dates))) {
     stop('`data$timeseries$date` must match `target_infection_dates` ',
-         'exactly for an `epiwave_greta_timeseries` object.')
+         'exactly for a `greta_timeseries` object.')
   }
 
   ihr <- data$ihr

R/compute_flat_prior_inits.R

@@ -0,0 +1,70 @@
+#' Compute initial values for the flat_prior infection model
+#'
+#' @description GAM-smoothed initial values and the observable-date window
+#'  are only used by `infection_model_type = 'flat_prior'` -- GP-based
+#'  infection models (`gp_infections`/`gp_growth_rate`/
+#'  `gp_growth_rate_deriv`) never reference `observable_idx`/`inits_values`.
+#'  This is deliberately not computed during data preparation
+#'  (`prepare_observation_data()`/`define_observation_model()`/
+#'  `stack_jurisdictions()`), since doing so unconditionally would run an
+#'  `mgcv::gam()` fit per stream per jurisdiction for no purpose whenever a
+#'  GP model is used instead. `fit_waves()` calls this lazily, only inside
+#'  its `flat_prior` branch, on the fully stacked multi-jurisdiction object.
+#'
+#' @param observations a stacked observations object, as produced by
+#'  `stack_jurisdictions()`
+#'
+#' @return list with `incidence_observable` (`date x jurisdiction` logical
+#'  matrix, `TRUE` where at least one stream has nearby information to infer
+#'  an infection) and `incidence_observable_inits` (`date x jurisdiction`
+#'  numeric matrix, mean GAM-smoothed initial value across streams)
+#' @noRd
+compute_flat_prior_inits <- function (observations) {
+
+  target_infection_dates <- observations$target_infection_dates
+  jurisdictions <- observations$target_jurisdictions
+  n_dates <- length(target_infection_dates)
+  stream_ids <- names(observations$observation_model_data)
+
+  per_jurisdiction <- lapply(seq_along(jurisdictions), function (j) {
+
+    stream_inits <- lapply(stream_ids, function (stream_id) {
+      stream <- observations$observation_model_data[[stream_id]]
+      # prop_mat_raw, not prop_mat: inits must use the pre-DOW-correction
+      # proportion (see stack_stream()) -- using the corrected one would make
+      # this depend on a prior-predictive draw of the DOW effect instead of
+      # being a deterministic approximation
+      inits_by_jurisdiction(
+        stream$case_mat[, j],
+        stream$delays[[j]],
+        stream$prop_mat_raw[, j],
+        target_infection_dates)
+    })
+
+    observable_idx <- rep(FALSE, n_dates)
+    inits_stream_mat <- matrix(NA_real_, n_dates, length(stream_ids))
+
+    for (s in seq_along(stream_inits)) {
+      idx <- stream_inits[[s]]$observable_idx
+      observable_idx[idx] <- TRUE
+      inits_stream_mat[idx, s] <- stream_inits[[s]]$inits_values
+    }
+
+    list(
+      observable_idx = observable_idx,
+      inits_values = rowMeans(inits_stream_mat, na.rm = TRUE)
+    )
+  })
+
+  incidence_observable <- do.call(
+    cbind, lapply(per_jurisdiction, `[[`, "observable_idx"))
+  incidence_observable_inits <- do.call(
+    cbind, lapply(per_jurisdiction, `[[`, "inits_values"))
+  colnames(incidence_observable) <- jurisdictions
+  colnames(incidence_observable_inits) <- jurisdictions
+
+  list(
+    incidence_observable = incidence_observable,
+    incidence_observable_inits = incidence_observable_inits
+  )
+}

R/create_epiwave_timeseries.R

@@ -1,51 +1,3 @@
-#' Expand constant value into a long tibble
-#'
-#' @description The epiwave model functions expect data in a long format,
-#'   structured to have a value for every date. This function creates a
-#'   tibble of this structure from a single value that is replicated in each
-#'   cell.
-#'
-#' @param dates infection dates sequence
-#' @param value value to be replicated in each cell
-#'
-#' @importFrom tibble tibble
-#'
-#' @return a long tibble with the constant value replicated across all
-#'   dates
-#' @export
-create_epiwave_timeseries <- function(dates,
-                                      value) {
-  long_combined <- tibble::tibble(date = dates,
-                                  value = value)
-
-  class(long_combined) <- c("epiwave_timeseries", class(long_combined))
-  long_combined
-}
-
-#' Expand a fixed value into a long tibble
-#'
-#' @description The epiwave model functions expect data in a long format,
-#'   structured to have a value for every date. This function creates a
-#'   tibble of this structure from a single fixed numeric value that is
-#'   replicated in each cell.
-#'
-#' @param dates infection dates sequence
-#' @param value a fixed numeric value to be replicated in each cell
-#'
-#' @return a long tibble with the fixed value replicated across all
-#'   dates, with class `epiwave_fixed_timeseries`
-#' @export
-create_epiwave_fixed_timeseries <- function(dates,
-                                            value) {
-  timeseries <- create_epiwave_timeseries(
-    dates = dates,
-    value = value
-  )
-
-  class(timeseries) <- c("epiwave_fixed_timeseries", class(timeseries))
-  timeseries
-}
-
 #' Create a greta-compatible timeseries object
 #'
 #' @description The epiwave model functions expect data in a long format,
@@ -61,13 +13,13 @@ create_epiwave_fixed_timeseries <- function(dates,
 #'
 #' @importFrom tibble tibble
 #'
-#' @return a list with class `epiwave_greta_timeseries` containing
-#'   components `timeseries` (a tibble of dates) and `ihr` (a greta array of
-#'   implied hospitalisation rates)
+#' @return a list with class `greta_timeseries` containing components
+#'   `timeseries` (a tibble of dates) and `ihr` (a greta array of implied
+#'   hospitalisation rates)
 #' @export
-create_epiwave_greta_timeseries <- function(dates,
-                                            car,
-                                            chr_prior) {
+as_greta_timeseries <- function(dates,
+                                car,
+                                chr_prior) {
   long_unique <- tibble::tibble(date = dates)
 
   dim(chr_prior) <- length(dates)
@@ -81,8 +33,47 @@ create_epiwave_greta_timeseries <- function(dates,
   long_combined <- list(timeseries = long_unique,
                         ihr = ihr_greta)
 
-  class(long_combined) <- c("epiwave_greta_timeseries",
+  class(long_combined) <- c("greta_timeseries",
                             "epiwave_timeseries",
                             class(long_combined))
   long_combined
 }
+
+#' Coerce a date/value table to an epiwave_fixed_timeseries object
+#'
+#' @description Observation data (case counts, hospitalisation counts) is
+#'   naturally a table of `date`/`value` pairs with real per-date
+#'   observations, rather than a
+#'   single value replicated across dates -- unlike delay distributions or
+#'   proportions, which usually apply uniformly. This function takes such a
+#'   table (with genuinely partial date coverage, e.g. missing weekends, is
+#'   fine -- gaps are handled downstream by `as_matrix()`) and classes it as
+#'   an `epiwave_fixed_timeseries` object, so it doesn't need to be
+#'   hand-classed with `class(x) <- c(...)`.
+#'
+#'   Already-classed `epiwave_timeseries` objects and plain numeric
+#'   values/vectors are returned unchanged, so this is safe to call
+#'   unconditionally on anything accepted as observation data.
+#'
+#' @param data an `epiwave_timeseries` object, a numeric value/vector, or a
+#'   data.frame/tibble with `date` and `value` columns
+#'
+#' @return an `epiwave_timeseries` or numeric object, ready for `as_matrix()`
+#' @noRd
+as_epiwave_timeseries <- function(data) {
+
+  if (inherits(data, "epiwave_timeseries") || is.numeric(data)) {
+    return(data)
+  }
+
+  if (!all(c("date", "value") %in% names(data))) {
+    stop("`data` must have `date` and `value` columns to be used as ",
+         "observation data")
+  }
+
+  data <- tibble::as_tibble(data[c("date", "value")])
+  data$date <- as.Date(data$date)
+
+  class(data) <- c("epiwave_fixed_timeseries", "epiwave_timeseries", class(data))
+  data
+}

R/define_observation_data.R

@@ -5,10 +5,11 @@
 #'  multiple jurisdictions later via `stack_jurisdictions()`.
 #'
 #' @param timeseries_data timeseries data for data of interest, for one
-#'  jurisdiction
-#' @param delay_from_infection a `discrete_pmf`/`discrete_weights` object
-#'  (replicated across dates), or an already time-varying
-#'  `discrete_pmf_series`/`discrete_weights_series` object
+#'  jurisdiction. A plain data.frame/tibble with `date` and `value` columns
+#'  is fine (gaps in date coverage are fine too) -- it's coerced
+#'  automatically, no need to pre-class it
+#' @param delay_from_infection a `discrete_pmf` object (replicated across
+#'  dates), or an already time-varying `discrete_pmf_series` object
 #' @param proportion_infections proportion data
 #' @param dow_model logical indicating whether to apply a DOW
 #'
@@ -26,42 +27,3 @@ define_observation_data <- function (timeseries_data,
   out
 
 }
-
-#' Define seroprevalence observation data
-#'
-#' @description Bundle one jurisdiction's seroprevalence survey data. Call
-#'  this once per jurisdiction; combine multiple jurisdictions later via
-#'  `stack_jurisdictions()`.
-#'
-#' @param timeseries_data seroprevalence survey timeseries data, for one
-#'  jurisdiction
-#' @param total_pop population size of this jurisdiction (a single numeric
-#'  value)
-#' @param size_vec sample size of the seroprevalence survey, for one
-#'  jurisdiction
-#' @param delay_from_infection typically a `discrete_weights`/
-#'  `discrete_weights_series` object (e.g. probability of testing
-#'  seropositive by day since infection, which need not sum to 1 -- unlike
-#'  case/hospitalisation notification, seroconversion is usually persistent
-#'  rather than a one-time event). A `discrete_pmf`/`discrete_pmf_series` is
-#'  also accepted if a normalised delay is more appropriate.
-#' @param proportion_infections proportion data
-#' @param dow_model logical indicating whether to apply a DOW
-#'
-#' @return list of observation data for one data type
-define_sero_data <- function (timeseries_data,
-                              total_pop,
-                              size_vec,
-                              delay_from_infection,
-                              proportion_infections,
-                              dow_model = FALSE) {
-
-  out <- list(timeseries_data = timeseries_data,
-              total_pop = total_pop,
-              size_vec = size_vec,
-              delay_from_infection = delay_from_infection,
-              proportion_infections = proportion_infections,
-              dow_model = dow_model)
-  out
-
-}