R/estimate_R_cl_single.R

In ern: Effective Reproduction Number Estimation

#' @title Correct underreporting by scaling up
#'
#' @param reports.daily dataframe of daily reported cases. must at least include `value` column with counts.
#' @param reporting.fraction numeric. proportion of incidence that is reported.
#' @keywords internal 
correct_underreporting <- function(
    reports.daily,
    reporting.fraction
){
  (reports.daily
   |> dplyr::mutate(value = value/reporting.fraction)
  )
}

#' @title Infer incidence from reports via a series of deconvolutions
#'
#' @param reports.daily Data frame. Daily report counts. Includes at least `date` and `value` columns.
#' @param reporting.delay List. Parameters for a single reporting delay distribution, as generated by [sample_a_dist()].
#' @param incubation.period List. Parameters for a single incubation period distribution, as generated by [sample_a_dist()].
#' @inheritParams deconv
#' @inheritParams estimate_R_cl
#'
#' @keywords internal
reports_to_incidence <- function(
    reports.daily,
    reporting.delay,
    incubation.period,
    max.iter = 9,
    silent = FALSE
){

  # reports -> onsets
  # - - - - - - - - - - - - - - - - -

  # Deconvolving inferred daily reports with reporting delay
  # distribution to get daily onsets.
  if(!silent){
    message("Deconvolution reporting delays...")
  }

  onsets <- deconv(
    counts   = reports.daily$value,
    dist     = get_discrete_dist(reporting.delay),
    max.iter = max.iter,
    silent   = silent
  )

  # onsets -> incidence
  # - - - - - - - - - - - - - - - - -

  # Deconvolving daily onsets with incubation period
  # distribution to get daily incidence
  if(!silent){
    message("Deconvolution incubation period...")
  }
  incidence <- deconv(
    counts   = onsets$y,
    dist     = get_discrete_dist(incubation.period),
    max.iter = max.iter,
    silent   = silent
  )

  # - - - output

  date.lookup <- tibble::tibble(
    date = reports.daily$date,
    t = 1:nrow(reports.daily)
  )

  res = incidence |>
    # filter out first x days, where x is the sum of
    # the max reporting delay and the max incubation period
    # (i.e. disregard reports before a full observation period is complete)
    # |> dplyr::filter(t >= incubation.period$max + reporting.delay$max)
    dplyr::left_join(date.lookup, by = "t") |>
    dplyr::transmute(date, I = y) |>
    tibble::as_tibble()

  return(res)
}

#' @title Wrapper for deconvolution with a given distribution
#'
#' @param counts Numeric. Vector of daily counts.
#' @param dist Numeric. Vector of discrete distributions (daily, truncated) with which we're deconvoluting counts, _e.g._, produced by [`get_discrete_dist()`]
#' @param max.iter Numeric. Maximum number of Richardson-Lucy iterations.
#' @inheritParams estimate_R_cl
#'
#' @keywords internal
deconv <- function(
    counts,
    dist,
    max.iter = 10,
    silent = FALSE
){
  # check args
  # - - - - - - - - - - - - - - - - -
  check_for_deconv(counts, dist)

  # perform Richardson-Lucy deconvolution
  (deconvolution_RL(
    observed = counts,
    times = 1:length(counts),
    p_delay = dist,
    max_iter = max.iter,
    verbose = !silent
  )
  |> tidyr::drop_na()
  |> dplyr::rename(
    t = time,
    y = RL_result)
  )
}




#' @title A single realization of the Rt estimate
#'
#' @param cl.daily Dataframe of inferred daily incidence.
#' @inheritParams estimate_R_cl
#' @keywords internal
#'
estimate_R_cl_single <- function(
    cl.daily,
    dist.repfrac,
    dist.repdelay,
    dist.incub,
    dist.gi,
    prm.R,
    silent = FALSE,
    RL.max.iter = 10
){
  
  # sample one realization of reports.daily (smoothed)
  id.list <- unique(cl.daily$id)
  the_id  <- sample(id.list, size = 1)
  df.draw <- (cl.daily
              |> dplyr::filter(id == the_id)
  )
  
  # sample reporting fraction
  reporting.fraction.draw <- sample_from_dist(
    n = 1,
    params = dist.repfrac
  )
  
  # sample reporting delay distribution
  reporting.delay <- sample_a_dist(dist.repdelay)
  
  # sample incubation period
  incubation.period <- sample_a_dist(dist.incub)
  
  # sample generation interval
  generation.interval <- sample_a_dist(dist.gi)
  
  # correct for underreporting
  reports.daily.scaled <- correct_underreporting(
    df.draw,
    reporting.fraction.draw
  )
  
  # reports deconvoluted with reporting delay distribution
  # and then with incubation period distribution
  incidence <- reports_to_incidence(
    reports.daily = reports.daily.scaled,
    reporting.delay,
    incubation.period,
    silent = silent,
    max.iter = RL.max.iter  ) |>
    # attach time index to incidence
    (\(x){dplyr::mutate(x, t = 1:nrow(x))})()
  
  
  # The RL deconvolution can be unstable at
  # the start of the deconvoluted function,
  # so we remove the first points based on 
  # the mean of the kernels (this is heuristic):
  incidence = dplyr::filter(incidence,
                            t > (dist.repdelay$mean + dist.incub$mean))
  
  # Estimate Rt
  res = incidence_to_R(
    incidence,
    generation.interval,
    prm.R)
  
  return(res)
}