Add some model components

NAMESPACE

@@ -28,12 +28,8 @@ export(get_ww_data_indices)
 export(get_ww_data_sizes)
 export(get_ww_values)
 export(indicate_ww_exclusions)
-export(make_hospital_onset_delay_pmf)
-export(make_incubation_period_pmf)
-export(make_reporting_delay_pmf)
 export(preprocess_hosp_data)
 export(preprocess_ww_data)
-export(simulate_double_censored_pmf)
 export(to_simplex)
 export(validate_paramlist)
 export(wwinference)

R/figures.R

@@ -3,7 +3,9 @@
 #' @param draws A dataframe containing the posterior draws with the data joined
 #' to it. This is the `draws_df` output of a call to [wwinference()]
 #' @param count_data_eval A dataframe containing the count data we will
-#' evaluate the forecasts against.
+#' evaluate the forecasts against. Must contain the columns `date` and
+#' a column indicating the count data to evaluate against, with the name
+#' of that column specified as the `count_data_eval_col_name`
 #' @param count_data_eval_col_name string indicating the name of the count
 #' data to evaluate against the forecasted count data
 #' @param forecast_date A string indicating the date we made the forecast, for

R/utils.R

@@ -141,3 +141,33 @@ convert_to_logsd <- function(mean, sd) {
 to_simplex <- function(vector) {
   return(vector / sum(vector))
 }
+
+#' Drop the first element of a simplex and re-normalize the result to sum to 1.
+#'
+#' When this vector corresponds to the generation interval distribution, we
+#' want to drop this first bin. The renewal equation assumes that same-day
+#' infection and onward transmission does not occur, and we assume
+#' everything is 1 indexed not 0 indeced. We need to
+#' manually drop the first element from the PMF vector.
+#'
+#' @param x A numeric vector, sums to 1. Corresponds to a discretized PDF or PMF
+#'   (usually the GI distribution).
+#'
+#' @return A numeric vector, sums to 1.
+#' @export
+#' @examples
+#' pmf_orig <- c(0.1, 0.1, 0.1, 0.7)
+#' pmf_trunc <- drop_first_and_renormalize(pmf_orig)
+drop_first_and_renormalize <- function(x) {
+  # Check input sums to 1
+  stopifnot(abs(sum(x) - 1) < 1e-8)
+  # Drop and renormalize
+  y <- x[2:length(x)] / sum(x[2:length(x)])
+  vec_outside_tol <- abs(sum(y) - 1L) > 1e-10
+  # Normalize until within tolerance
+  while (vec_outside_tol) {
+    y <- y / sum(y)
+    vec_outside_tol <- abs(sum(y) - 1L) > 1e-10
+  }
+  return(y)
+}

data-raw/covid_pmfs.R

@@ -1,9 +1,116 @@
+# Load in the parameters
 params <- get_params(
   system.file("extdata", "example_params.toml",
     package = "wwinference"
   )
 )
 
+
+
+
+# Simulate daily double censored PMF. From {epinowcast}:
+# https://package.epinowcast.org/dev/reference/simulate_double_censored_pmf.html #nolint
+#
+# This function simulates the probability mass function of a  daily
+# double-censored process. The process involves two distributions: a primary
+# distribution which represents the censoring process for the primary event
+# and another distribution (which is offset by the primary).
+#
+# Based off of:
+# https://www.medrxiv.org/content/10.1101/2024.01.12.24301247v1
+simulate_double_censored_pmf <- function(
+    max, fun_primary = stats::runif, primary_args = list(),
+    fun_dist = stats::rlnorm,
+    dist_args = list(...), n = 1e6, ...) {
+  primary <- do.call(fun_primary, c(list(n), primary_args))
+  secondary <- primary + do.call(fun_dist, c(list(n), dist_args))
+  delay <- floor(secondary) - floor(primary)
+  if (missing(max)) {
+    max <- base::max(delay)
+  }
+  cdf <- ecdf(delay)(0:max)
+  pmf <- c(cdf[1], diff(cdf))
+  vec_outside_tol <- abs(sum(pmf) - 1L) > 1e-10
+  while (vec_outside_tol) {
+    pmf <- pmf / sum(pmf)
+    vec_outside_tol <- abs(sum(pmf) - 1L) > 1e-10
+  }
+  return(pmf)
+}
+
+
+# This retruns an incubation period pmf corresponding to the incubation period
+# for COVID after Omicron used in
+# Park et al 2023. These estimates are from early Omicron.
+
+make_incubation_period_pmf <- function(backward_scale = 3.60,
+                                       backward_shape = 1.50,
+                                       r = 0.15) {
+  # From: Park, Sang Woo, et al. "Inferring the differences in incubation-period
+  # and generation-interval distributions of the Delta and Omicron variants of
+  # SARS-CoV-2." Proceedings of the National Academy of Sciences 120.22 (2023):
+  # e2221887120.
+
+  # "However, when we account for growth-rate differences and reestimate the
+  # forward incubation periods, we find that both variants have similar
+  # incubation-period distributions with a mean of 4.1 d (95% CI: 3.8 to 4.4 d)
+  # for the Delta variant and 4.2 d (95% CI: 3.6 to 4.9 d) for the Omicron
+  # variant Fig. 3B)."
+
+  # Fits a Weibull to the data
+
+  # Relies on fundamental assumption about epidemic growth rate.
+
+
+  discr_gr_adj_weibull <- tibble::tibble(
+    time = seq(0, 23, by = 1), # 23 seems to get most of the distribution mass
+    density0 = dweibull(time,
+      shape = backward_shape,
+      scale = backward_scale
+    ) * exp(r * time)
+  )
+
+  inc_period_pmf <- to_simplex(discr_gr_adj_weibull$density0)
+  return(inc_period_pmf)
+}
+
+# Makes the hospital onset delay pmf using the parameter estimates based on
+# Danache et al linelist data from symptom onset to hospital
+# admission. See below:
+# https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0261428
+
+make_hospital_onset_delay_pmf <- function(neg_binom_mu = 6.98665,
+                                          neg_binom_size = 2.490848) {
+  density <- dnbinom(
+    x = seq(0, 30, 1),
+    mu = neg_binom_mu, size = neg_binom_size
+  )
+  hosp_onset_delay_pmf <- density / sum(density)
+
+  return(hosp_onset_delay_pmf)
+}
+
+
+# Makes the reporting delay pmf by convolving the incubation period pmf with
+# the symptom to hospital admission pmf and normalizing
+
+make_reporting_delay_pmf <- function(incubation_period_pmf,
+                                     hospital_onset_delay_pmf) {
+  pmfs <- list(
+    "incubation_period" = incubation_period_pmf,
+    "hosp_onset_delay" = hospital_onset_delay_pmf
+  )
+
+  infection_to_hosp_delay_pmf <- add_pmfs(pmfs) |>
+    to_simplex()
+  return(infection_to_hosp_delay_pmf)
+}
+
+
+
+
+
+# Put it all together
 generation_interval <- withr::with_seed(42, {
   wwinference::simulate_double_censored_pmf(
     max = params$gt_max, meanlog = params$mu_gi,