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