Make target_infection_dates emergent instead of a required input

.gitignore

@@ -10,3 +10,4 @@ figures/
 simdata/
 .positai
 README.html
+Rplots.pdf

NAMESPACE

@@ -1,8 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(as_matrix,epiwave_timeseries)
-S3method(as_matrix,greta_timeseries)
 S3method(as_matrix,numeric)
+S3method(print,epiwave_observation_model)
+S3method(print,epiwave_stacked_observations)
 export(as_greta_timeseries)
 export(as_matrix)
 export(compute_reff)
@@ -16,15 +17,20 @@ export(fit_waves)
 export(implement_day_of_week)
 export(new_convolution_matrix)
 export(plot_infection_traj)
+export(plot_observation_coverage)
 export(stack_jurisdictions)
 importFrom(cowplot,panel_border)
 importFrom(cowplot,theme_cowplot)
 importFrom(dplyr,bind_rows)
 importFrom(dplyr,filter)
 importFrom(dplyr,mutate)
+importFrom(ggplot2,aes)
 importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_line)
+importFrom(ggplot2,geom_tile)
 importFrom(ggplot2,ggplot)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,scale_fill_manual)
 importFrom(ggplot2,scale_x_date)
 importFrom(ggplot2,theme)
 importFrom(greta,"%*%")
@@ -40,7 +46,6 @@ importFrom(greta.gp,gp)
 importFrom(greta.gp,mat52)
 importFrom(mgcv,gam)
 importFrom(mgcv,predict.gam)
-importFrom(tibble,tibble)
 importFrom(tidyr,any_of)
 importFrom(tidyr,pivot_longer)
 importFrom(tidyr,starts_with)

R/as_matrix.R

@@ -49,17 +49,3 @@ as_matrix.epiwave_timeseries <- function (data, target_infection_dates, ...) {
   names(out) <- as.character(target_dates)
   out
 }
-
-#' @export
-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 a `greta_timeseries` object.')
-  }
-
-  ihr <- data$ihr
-  dim(ihr) <- c(length(target_infection_dates), 1)
-
-  ihr
-}

R/create_epiwave_timeseries.R

@@ -1,42 +1,26 @@
-#' Create a greta-compatible timeseries object
+#' Bundle a case ascertainment rate and case-hospitalisation-rate prior
 #'
-#' @description The epiwave model functions expect data in a long format,
-#'   structured to have a value for every date. This function creates a
-#'   greta-compatible timeseries from a case ascertainment rate and a
-#'   hospitalisation rate prior.
+#' @description Builds a greta-backed `proportion_infections` input (e.g.
+#'   the implied hospitalisation rate, IHR = CAR x CHR) without needing a
+#'   date axis up front. The greta array this eventually becomes has to be
+#'   dimensioned to `target_infection_dates`, but that axis is emergent --
+#'   derived later, from the data, by `stack_jurisdictions()` -- so the
+#'   actual `car * chr_prior` multiplication and dimensioning is deferred
+#'   until then, once the axis is known.
 #'
-#' @param dates infection dates sequence
-#' @param car case ascertainment rate; either a numeric value or an
-#'   `epiwave_timeseries` object
+#' @param car case ascertainment rate; a fixed numeric value
 #' @param chr_prior a greta array representing the prior distribution for
 #'   the case hospitalisation rate
 #'
-#' @importFrom tibble tibble
-#'
-#' @return a list with class `greta_timeseries` containing components
-#'   `timeseries` (a tibble of dates) and `ihr` (a greta array of implied
-#'   hospitalisation rates)
+#' @return a list with class `greta_proportion` containing `car` and
+#'   `chr_prior`, unresolved until `stack_jurisdictions()`
 #' @export
-as_greta_timeseries <- function(dates,
-                                car,
+as_greta_timeseries <- function(car,
                                 chr_prior) {
-  long_unique <- tibble::tibble(date = dates)
-
-  dim(chr_prior) <- length(dates)
-
-  if ("epiwave_timeseries" %in% class(car)) {
-    car <- car$value
-  }
-
-  ihr_greta <- car * chr_prior
-
-  long_combined <- list(timeseries = long_unique,
-                        ihr = ihr_greta)
 
-  class(long_combined) <- c("greta_timeseries",
-                            "epiwave_timeseries",
-                            class(long_combined))
-  long_combined
+  out <- list(car = car, chr_prior = chr_prior)
+  class(out) <- c("greta_proportion", class(out))
+  out
 }
 
 #' Coerce a date/value table to an epiwave_fixed_timeseries object

R/define_observation_model.R

@@ -7,24 +7,33 @@
 #' combine multiple jurisdictions, pass several of these to
 #' `stack_jurisdictions()` first.
 #'
-#' @param target_infection_dates sequence of infection dates
+#' No `target_infection_dates` argument is needed here -- the date axis is
+#' emergent, derived later (by `stack_jurisdictions()`) from every stream's
+#' own data and delay distributions. This function just aggregates each
+#' stream's implied date range into one range for the whole jurisdiction.
+#'
 #' @param ... observation data sets for this jurisdiction (as returned by
 #'  `define_observation_data()`), named by stream
 #'
 #' @return list describing one jurisdiction's observation model, with class
 #'  `epiwave_observation_model`
 #' @export
 #'
-define_observation_model <- function (target_infection_dates = NULL, ...) {
+define_observation_model <- function (...) {
 
   observation_list <- list(...)
 
   prepared_observation_model_data <- lapply(observation_list,
-                                            prepare_observation_data,
-                                            target_infection_dates)
+                                            prepare_observation_data)
+
+  implied_ranges <- lapply(prepared_observation_model_data,
+                           function(x) x$implied_range)
+  starts <- do.call(c, lapply(implied_ranges, `[`, 1))
+  ends <- do.call(c, lapply(implied_ranges, `[`, 2))
+  implied_range <- c(min(starts), max(ends))
 
   out <- list(observation_model_data = prepared_observation_model_data,
-             target_infection_dates = target_infection_dates)
+             implied_range = implied_range)
 
   class(out) <- c("epiwave_observation_model", class(out))
   out

R/epiwave-package.R

@@ -11,6 +11,8 @@ if (getRversion() >= "2.15.1") {
 globalVariables(
   c(
     "x",
-    "z"
+    "z",
+    "jurisdiction",
+    "observed"
   )
 )

R/fit_waves.R

@@ -30,7 +30,13 @@
 #' @param n_samples number of samples after warmup
 #' @param n_extra_samples number of extra samples if original run didn't converge
 #'
-#' @return list of infection_model, fit, infection_days, and jurisdictions
+#' @return list of infection_model, fit, infection_days, jurisdictions,
+#'   observation_models, incidence_greta_arrays, and observation_model_data
+#'   (the resolved, axis-aligned per-stream data -- e.g. `prop_mat` for a
+#'   stream with a greta-backed `proportion_infections` like IHR-from-CHR is
+#'   only ever computed once the emergent axis is known, so it's exposed
+#'   here rather than remaining accessible on the object originally passed
+#'   to `as_greta_timeseries()`)
 #' @export
 #'
 fit_waves <- function (observations,
@@ -123,7 +129,8 @@ fit_waves <- function (observations,
     infection_days = target_infection_dates,
     jurisdictions = jurisdictions,
     observation_models = observation_models,
-    incidence_greta_arrays = incidence_greta_arrays)
+    incidence_greta_arrays = incidence_greta_arrays,
+    observation_model_data = observation_model_data)
 
   return(fit_output)
 

R/plot_observation_coverage.R

@@ -0,0 +1,122 @@
+#' Visualise observation data coverage against the derived date axis
+#'
+#' @description Since `target_infection_dates` is emergent rather than
+#'  user-supplied (see `stack_jurisdictions()`), it's worth being able to
+#'  see what actually got derived and how each stream/jurisdiction's data
+#'  aligns to it before committing to a (potentially slow) `fit_waves()`
+#'  run -- this is the visual form of the alignment checks used throughout
+#'  the package's own tests.
+#'
+#' @param observations either an already-stacked `epiwave_stacked_observations`
+#'  object (output of `stack_jurisdictions()`), or a single jurisdiction's raw
+#'  `epiwave_observation_model` (output of `define_observation_model()`),
+#'  which is stacked on its own to derive its axis first
+#'
+#' @importFrom dplyr bind_rows
+#' @importFrom tidyr pivot_longer
+#' @importFrom ggplot2 ggplot aes geom_tile facet_wrap scale_fill_manual
+#'  scale_x_date labs
+#' @importFrom cowplot theme_cowplot panel_border
+#'
+#' @return a ggplot object: a date x jurisdiction coverage tile per stream
+#' @export
+plot_observation_coverage <- function (observations) {
+
+  if (!inherits(observations, 'epiwave_stacked_observations')) {
+    observations <- stack_jurisdictions_list(list(observations))
+  }
+
+  target_infection_dates <- observations$target_infection_dates
+
+  coverage_df <- dplyr::bind_rows(lapply(
+    names(observations$observation_model_data),
+    function (stream_id) {
+      case_mat <- observations$observation_model_data[[stream_id]]$case_mat
+      # go via column names (jurisdiction labels), not positional/vector
+      # order, to avoid silently mismatching dates to jurisdictions
+      observed_df <- as.data.frame(!is.na(case_mat), check.names = FALSE)
+      observed_df$date <- target_infection_dates
+      observed_long <- tidyr::pivot_longer(
+        observed_df,
+        cols = -date,
+        names_to = "jurisdiction",
+        values_to = "observed")
+      observed_long$stream <- stream_id
+      observed_long
+    }))
+
+  ggplot2::ggplot(
+    coverage_df,
+    ggplot2::aes(x = date, y = jurisdiction, fill = observed)) +
+    ggplot2::geom_tile() +
+    ggplot2::facet_wrap(~stream, ncol = 1) +
+    ggplot2::scale_fill_manual(values = c(`TRUE` = "steelblue", `FALSE` = "grey90")) +
+    cowplot::theme_cowplot() +
+    cowplot::panel_border(remove = TRUE) +
+    ggplot2::scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
+    ggplot2::labs(x = NULL, y = NULL, fill = "observed")
+}
+
+#' Print a summary of a stacked observation model
+#'
+#' @description Reports the derived `target_infection_dates` axis and, per
+#'  stream and jurisdiction, the observed date range and coverage against
+#'  that axis -- a quick text-only check of what `stack_jurisdictions()`
+#'  actually derived, for when a full `plot_observation_coverage()` isn't
+#'  needed.
+#'
+#' @param x an `epiwave_stacked_observations` object
+#' @param ... unused, present for consistency with the `print()` generic
+#'
+#' @return `x`, invisibly
+#' @export
+print.epiwave_stacked_observations <- function (x, ...) {
+
+  target_infection_dates <- x$target_infection_dates
+
+  cat("<epiwave_stacked_observations>\n")
+  cat(sprintf(
+    "target_infection_dates: %s to %s (%d days)\n",
+    format(min(target_infection_dates)),
+    format(max(target_infection_dates)),
+    length(target_infection_dates)))
+  cat("jurisdictions:", paste(x$target_jurisdictions, collapse = ", "), "\n")
+
+  for (stream_id in names(x$observation_model_data)) {
+    cat("\n", stream_id, ":\n", sep = "")
+    case_mat <- x$observation_model_data[[stream_id]]$case_mat
+    for (jurisdiction in x$target_jurisdictions) {
+      observed <- !is.na(case_mat[, jurisdiction])
+      if (any(observed)) {
+        observed_dates <- target_infection_dates[observed]
+        cat(sprintf(
+          "  %s: %s to %s (%d of %d days observed)\n",
+          jurisdiction,
+          format(min(observed_dates)),
+          format(max(observed_dates)),
+          sum(observed),
+          length(target_infection_dates)))
+      } else {
+        cat(sprintf("  %s: no observed data\n", jurisdiction))
+      }
+    }
+  }
+
+  invisible(x)
+}
+
+#' Print a summary of a single jurisdiction's observation model
+#'
+#' @description Derives this jurisdiction's own axis (as `fit_waves()` would
+#'  for a single-jurisdiction fit) and reports it the same way as
+#'  `print.epiwave_stacked_observations()`.
+#'
+#' @param x an `epiwave_observation_model` object
+#' @param ... unused, present for consistency with the `print()` generic
+#'
+#' @return `x`, invisibly
+#' @export
+print.epiwave_observation_model <- function (x, ...) {
+  print(stack_jurisdictions_list(list(x)))
+  invisible(x)
+}