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