nmode's Git Repositories - Rnaught/blob - man/WP.Rd

   1 % Generated by roxygen2: do not edit by hand
   2 % Please edit documentation in R/WP.R
   3 \name{WP}
   4 \alias{WP}
   5 \title{WP method}
   6 \usage{
   7 WP(
   8   NT,
   9   mu = NA,
  10   search = list(B = 100, shape.max = 10, scale.max = 10),
  11   tol = 0.999
  12 )
  13 }
  14 \arguments{
  15 \item{NT}{Vector of case counts.}
  16
  17 \item{mu}{Mean of the serial distribution (needs to match case counts in time units; for example, if case
  18 counts are weekly and the serial distribution has a mean of seven days, then \code{mu} should be
  19 set to one). The default value of \code{mu} is set to \code{NA}.}
  20
  21 \item{search}{List of default values for the grid search algorithm. The list includes three elements: the
  22 first is \code{B}, which is the length of the grid in one dimension; the second is
  23 \code{scale.max}, which is the largest possible value of the scale parameter; and the third
  24 is \code{shape.max}, which is the largest possible value of the shape parameter. Defaults to
  25 \code{B=100, scale.max=10, shape.max=10}. For both shape and scale, the smallest possible
  26 value is 1/\code{B}.}
  27
  28 \item{tol}{Cutoff value for cumulative distribution function of the pre-discretization gamma serial
  29 distribution. Defaults to 0.999 (i.e. in the discretization, the maximum is chosen such that the
  30 original gamma distribution has cumulative probability of no more than 0.999 at this maximum).}
  31 }
  32 \value{
  33 \code{WP} returns a list containing the following components:  \code{Rhat} is the estimate of R0,
  34         and \code{SD} is either the discretized serial distribution (if \code{mu} is not \code{NA}), or the
  35         estimated discretized serial distribution (if \code{mu} is \code{NA}). The list also returns the
  36         variable \code{check}, which is equal to the number of non-unique maximum likelihood estimators.
  37         The serial distribution \code{SD} is returned as a list made up of \code{supp} (the support of
  38         the distribution) and \code{pmf} (the probability mass function).
  39 }
  40 \description{
  41 This function implements an R0 estimation due to White and Pagano (Statistics in Medicine, 2008).
  42 The method is based on maximum likelihood estimation in a Poisson transmission model.
  43 See details for important implementation notes.
  44 }
  45 \details{
  46 This method is based on a Poisson transmission model, and hence may be most most valid at the beginning
  47 of an epidemic. In their model, the serial distribution is assumed to be discrete with a finite number
  48 of posible values. In this implementation, if \code{mu} is not {NA}, the serial distribution is taken to
  49 be a discretized version of a gamma distribution with mean \code{mu}, shape parameter one, and largest
  50 possible value based on parameter \code{tol}. When \code{mu} is \code{NA}, the function implements a
  51 grid search algorithm to find the maximum likelihood estimator over all possible gamma distributions
  52 with unknown mean and variance, restricting these to a prespecified grid (see \code{search} parameter).
  53
  54 When the serial distribution is known (i.e., \code{mu} is not \code{NA}), sensitivity testing of \code{mu}
  55 is strongly recommended. If the serial distribution is unknown (i.e., \code{mu} is \code{NA}), the
  56 likelihood function can be flat near the maximum, resulting in numerical instability of the optimizer.
  57 When \code{mu} is \code{NA}, the implementation takes considerably longer to run. Users should be careful
  58 about units of time (e.g., are counts observed daily or weekly?) when implementing.
  59
  60 The model developed in White and Pagano (2008) is discrete, and hence the serial distribution is finite
  61 discrete. In our implementation, the input value \code{mu} is that of a continuous distribution. The
  62 algorithm discretizes this input when \code{mu} is not \code{NA}, and hence the mean of the serial
  63 distribution returned in the list \code{SD} will differ from \code{mu} somewhat. That is to say, if the
  64 user notices that the input \code{mu} and output mean of \code{SD} are different, this is to be expected,
  65 and is caused by the discretization.
  66 }
  67 \examples{
  68 ## ===================================================== ##
  69 ## Illustrate on weekly data                             ##
  70 ## ===================================================== ##
  71
  72 NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  73 ## obtain Rhat when serial distribution has mean of five days
  74 res1 <- WP(NT=NT, mu=5/7)
  75 res1$Rhat
  76 ## obtain Rhat when serial distribution has mean of three days
  77 res2    <- WP(NT=NT, mu=3/7)
  78 res2$Rhat
  79 ## obtain Rhat when serial distribution is unknown
  80 ## NOTE:  this implementation will take longer to run
  81 res3    <- WP(NT=NT)
  82 res3$Rhat
  83 ## find mean of estimated serial distribution
  84 serial <- res3$SD
  85 sum(serial$supp * serial$pmf)
  86
  87 ## ========================================================= ##
  88 ## Compute Rhat using only the first five weeks of data      ##
  89 ## ========================================================= ##
  90
  91 res4 <- WP(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days
  92 res4$Rhat
  93
  94 }