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

   1 % Generated by roxygen2: do not edit by hand
   2 % Please edit documentation in R/seqB.R
   3 \name{seqB}
   4 \alias{seqB}
   5 \title{seqB method}
   6 \usage{
   7 seqB(NT, mu, kappa = 20)
   8 }
   9 \arguments{
  10 \item{NT}{Vector of case counts.}
  11
  12 \item{mu}{Mean of the serial distribution. This needs to match case counts in
  13 time units. For example, if case counts are weekly and the serial
  14 distribution has a mean of seven days, then \code{mu} should be set
  15 to one. If case counts are daily and the serial distribution has a
  16 mean of seven days, then \code{mu} should be set to seven.}
  17
  18 \item{kappa}{Largest possible value of uniform prior (defaults to 20). This
  19 describes the prior belief on ranges of R0, and should be set to
  20 a higher value if R0 is believed to be larger.}
  21 }
  22 \value{
  23 \code{seqB} returns a list containing the following components:
  24         \code{Rhat} is the estimate of R0 (the posterior mean),
  25         \code{posterior} is the posterior distribution of R0 from which
  26         alternate estimates can be obtained (see examples), and \code{group}
  27         is an indicator variable (if \code{group == TRUE}, zero values of NT
  28         were input and grouping was done to obtain \code{Rhat}). The variable
  29         \code{posterior} is returned as a list made up of \code{supp} (the
  30         support of the distribution) and \code{pmf} (the probability mass
  31         function).
  32 }
  33 \description{
  34 This function implements a sequential Bayesian estimation method of R0 due to
  35 Bettencourt and Riberio (PloS One, 2008). See details for important
  36 implementation notes.
  37 }
  38 \details{
  39 The method sets a uniform prior distribution on R0 with possible values
  40 between zero and \code{kappa}, discretized to a fine grid. The distribution
  41 of R0 is then updated sequentially, with one update for each new case count
  42 observation. The final estimate of R0 is \code{Rhat}, the mean of the (last)
  43 posterior distribution. The prior distribution is the initial belief of the
  44 distribution of R0, which is the uninformative uniform distribution with
  45 values between zero and \code{kappa}. Users can change the value of
  46 \code{kappa} only (i.e., the prior distribution cannot be changed from the
  47 uniform). As more case counts are observed, the influence of the prior
  48 distribution should lessen on the final estimate \code{Rhat}.
  49
  50 This method is based on an approximation of the SIR model, which is most
  51 valid at the beginning of an epidemic. The method assumes that the mean of
  52 the serial distribution (sometimes called the serial interval) is known. The
  53 final estimate can be quite sensitive to this value, so sensitivity testing
  54 is strongly recommended. Users should be careful about units of time (e.g.,
  55 are counts observed daily or weekly?) when implementing.
  56
  57 Our code has been modified to provide an estimate even if case counts equal
  58 to zero are present in some time intervals. This is done by grouping the
  59 counts over such periods of time. Without grouping, and in the presence of
  60 zero counts, no estimate can be provided.
  61 }
  62 \examples{
  63 # Weekly data.
  64 NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  65
  66 ## Obtain R0 when the serial distribution has a mean of five days.
  67 res1 <- seqB(NT, mu = 5 / 7)
  68 res1$Rhat
  69
  70 ## Obtain R0 when the serial distribution has a mean of three days.
  71 res2 <- seqB(NT, mu = 3 / 7)
  72 res2$Rhat
  73
  74 # Compute posterior mode instead of posterior mean and plot.
  75
  76 Rpost <- res1$posterior
  77 loc <- which(Rpost$pmf == max(Rpost$pmf))
  78 Rpost$supp[loc] # Posterior mode.
  79 res1$Rhat # Compare with the posterior mean.
  80
  81 par(mfrow = c(2, 1), mar = c(2, 2, 1, 1))
  82
  83 plot(Rpost$supp, Rpost$pmf, col = "black", type = "l", xlab = "", ylab = "")
  84 abline(h = 1 / (20 / 0.01 + 1), col = "red")
  85 abline(v = res1$Rhat, col = "blue")
  86 abline(v = Rpost$supp[loc], col = "purple")
  87 legend("topright",
  88   legend = c("Prior", "Posterior", "Posterior mean", "Posterior mode"),
  89   col = c("red", "black", "blue", "purple"), lty = 1)
  90
  91 }