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. This needs to match case counts in time units. For example, if case counts
  17 #'           are 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 single value, the estimate of R0.
  21 #'
  22 #' @examples
  23 #' ## ===================================================== ##
  24 #' ## Illustrate on weekly data                             ##
  25 #' ## ===================================================== ##
  26 #'
  27 #' NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  28 #' ## obtain Rhat when serial distribution has mean of five days
  29 #' IDEA(NT=NT, mu=5/7)
  30 #' ## obtain Rhat when serial distribution has mean of three days
  31 #' IDEA(NT=NT, mu=3/7)
  32 #'
  33 #' ## ========================================================= ##
  34 #' ## Compute Rhat using only the first five weeks of data      ##
  35 #' ## ========================================================= ##
  36 #'
  37 #' IDEA(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days
  38 #'
  39 #' @export
  40 IDEA <- function(NT, mu) {
  41     if (length(NT) < 2)
  42         print("Warning: length of NT should be at least two.")
  43     else {
  44         NT <- as.numeric(NT)
  45         TT <- length(NT)
  46         s <- (1:TT) / mu
  47
  48         y1 <- log(NT) / s
  49         y2 <- s^2
  50         y3 <- log(NT)
  51
  52         IDEA1 <- sum(y2) * sum(y1) - sum(s) * sum(y3)
  53         IDEA2 <- TT * sum(y2) - sum(s)^2
  54         IDEA <- exp(IDEA1 / IDEA2)
  55
  56         return(IDEA)
  57     }
  58 }