nmode's Git Repositories - Rnaught/blob - R/IDEA.R

   1 #' IDEA method
   2 #'
   3 #' This function implements a least squares estimation method of R0 due to Fisman et al. (PloS One, 2013).
   4 #' See details for implementation notes.
   5 #'
   6 #' This method is closely related to that implemented in \code{ID}. The method is based on an incidence decay model.
   7 #' The estimate of R0 is the value which minimizes the sum of squares between observed case counts and cases counts
   8 #' expected under the model.
   9 #'
  10 #' This method is based on an approximation of the SIR model, which is most valid at the beginning of an epidemic.
  11 #' The method assumes that the mean of the serial distribution (sometimes called the serial interval) is known.
  12 #' The final estimate can be quite sensitive to this value, so sensitivity testing is strongly recommended.
  13 #' Users should be careful about units of time (e.g. are counts observed daily or weekly?) when implementing.
  14 #'
  15 #' @param NT Vector of case counts
  16 #' @param mu Mean of the serial distribution (needs to match case counts in time units; for example, if case counts are
  17 #'           weekly and the serial distribution has a mean of seven days, then \code{mu} should be set to one, if case
  18 #'           counts are daily and the serial distribution has a mean of seven days, then \code{mu} should be set to seven)
  19 #'
  20 #' @return \code{IDEA} returns a list containing the following components: \code{Rhat} is the estimate of R0 and
  21 #'         \code{inputs} is a list of the original input variables \code{NT, mu}.
  22 #'
  23 #' @examples
  24 #' ## ===================================================== ##
  25 #' ## Illustrate on weekly data                             ##
  26 #' ## ===================================================== ##
  27 #'
  28 #' NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  29 #' ## obtain Rhat when serial distribution has mean of five days
  30 #' IDEA(NT=NT, mu=5/7)
  31 #' ## obtain Rhat when serial distribution has mean of three days
  32 #' IDEA(NT=NT, mu=3/7)
  33 #'
  34 #' ## ========================================================= ##
  35 #' ## Compute Rhat using only the first five weeks of data      ##
  36 #' ## ========================================================= ##
  37 #'
  38 #' IDEA(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days
  39 #'
  40 #' @export
  41 IDEA <- function(NT, mu) {
  42     if (length(NT) < 2)
  43         print("Warning: length of NT should be at least two.")
  44     else {
  45         NT <- as.numeric(NT)
  46         TT <- length(NT)
  47         s <- (1:TT) / mu
  48
  49         y1 <- log(NT) / s
  50         y2 <- s^2
  51         y3 <- log(NT)
  52
  53         IDEA1 <- sum(y2) * sum(y1) - sum(s) * sum(y3)
  54         IDEA2 <- TT * sum(y2) - sum(s)^2
  55         IDEA <- exp(IDEA1 / IDEA2)
  56
  57         return(IDEA)
  58     }
  59 }