R/compute_force_infection.R

In reconhub/projections: Project Future Case Incidence

#' Calculate force of infection
#'
#' This function calculates the force of infection at time t+1, generated by
#' cases with onset at time 1, 2, ..., t, with reproduction numbers R_1, R_2,
#' ..., R_t. Calculations use matrices where independent simulations are stored
#' as separate columns.
#'
#' @param w a `numeric` containing numbers representing the PMF of the serial
#'   interval, starting at day 1, i.e. one day after symptom onset; the vector
#'   should always be at least as long as the largest value of `t` allowed
#'
#' @param cases a `matrix` of `integers` with `t` rows (each row is a date) and
#'   `n_sim` columns (each column is a separate simulation) indicating number of
#'   cases on a given day, in a given simulation
#'
#' @param R a `matrix` of `numeric` with `t` rows (each row is a date) and
#'   `n_sim` columns (each column is a separate simulation) indicating the
#'   reproduction number on a given day, in a given simulation
#'
#' @param t an `integer` indicating the simulation step: incidence will then be
#'   computed for `t+1` taking into account past cases and R from time point `1`
#'   until `t`
#'
#' @param instantaneous_R a boolean specifying whether to assume `R` is the case
#'   reproduction number (`instantaneous_R = FALSE`), or the
#'   instantaneous reproduction number (`instantaneous_R = TRUE`).
#'   If `instantaneous_R = FALSE` then values of `R` at time `t` will govern the
#'   mean number of secondary cases of all cases infected at time `t`,
#'   even if those secondary cases appear after `t`. In other words, `R`
#'   will characterise onwards transmission from infectors depending on their
#'   date of infection.
#'   If `instantaneous_R = TRUE` then values of `R` at time `t` will govern the
#'   mean number of secondary cases made at time `t` by all cases infected
#'   before `t`. In other words, `R` will characterise onwards transmission at
#'   a given time.
#'
#' @noRd

compute_force_infection <- function(w, cases, R, t, instantaneous_R) {
  rev_w <- rev(w)
  ws <- utils::tail(rev_w, t)

  cases <- cases[seq_len(t), , drop = FALSE]
  R <- R[seq_len(t), , drop = FALSE]

  if(!instantaneous_R) {
    lambda <- ws %*% (cases * R)
  } else {
    lambda <- (ws %*% cases) * R[t,]
  }
  as.vector(lambda)
}