nmode's Git Repositories - Rnaught/blob - vignettes/wp_serial.Rmd

   1 ---
   2 title: "White and Pagano: Utilizing the Serial Distribution"
   3 output: rmarkdown::html_vignette
   4 vignette: >
   5   %\VignetteIndexEntry{White and Pagano: Utilizing the Serial Distribution}
   6   %\VignetteEngine{knitr::rmarkdown}
   7   %\VignetteEncoding{UTF-8}
   8 ---
   9
  10 ```{r, include = FALSE}
  11 knitr::opts_chunk$set(
  12   collapse = TRUE,
  13   comment = "#>"
  14 )
  15 ```
  16
  17 ```{r setup, include = FALSE}
  18 library(Rnaught)
  19 ```
  20
  21 The serial distribution of an infectious disease is the distribution of the time
  22 from when an infectious individual -- the infector -- becomes symptomatic, to
  23 when another individual who is infected by the infector becomes symptomatic. The
  24 serial interval refers to a range of likely values from this distribution,
  25 although it is typically reported as the mean.
  26
  27 In the White and Pagano method, the serial distribution is assumed to be a
  28 discretized, finite version of a gamma distribution. Setting the parameter
  29 `serial` to `TRUE` causes this discretized distribution to be returned in
  30 addition to the estimate of R0. Furthermore, the method can be used whether or
  31 not the serial interval (specified as the parameter `mu`) is known. When `mu` is
  32 specified, it is taken to be the mean of a continuous gamma distribution (i.e.,
  33 before the discretization). As such, the mean computed from the returned serial
  34 distribution may differ slightly from `mu`:
  35
  36 ```{r}
  37 # Case counts.
  38 cases <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  39
  40 estimate <- wp(cases, mu = 3.333, serial = TRUE)
  41
  42 # `supp` is the support of the distribution, and `pmf` is its probability mass
  43 # function.
  44 sum(estimate$supp * estimate$pmf)
  45 ```
  46
  47 When `mu` is unspecified (left to its default value of `NA`), the method
  48 performs a maximum likelihood estimation over all (discretized) gamma
  49 distributions via a grid search, whose range of parameters are specified via
  50 `grid_length`, `max_shape` and `max_scale` (see `?wp` for more details). It is
  51 useful to return the estimated serial distribution in this case, as it can
  52 provide estimates of the serial interval when it is unknown:
  53
  54 ```{r}
  55 # The grid search parameters specified below are the default values.
  56 estimate <- wp(cases, serial = TRUE,
  57   grid_length = 100, max_shape = 10, max_scale = 10
  58 )
  59
  60 serial_mean <- sum(estimate$supp * estimate$pmf)
  61 serial_mean
  62
  63 # Compute the (discrete) median for an alternative estimate of the serial
  64 # interval.
  65 cdf <- cumsum(estimate$pmf)
  66 serial_med <- estimate$supp[which(cdf >= 0.5 & estimate$pmf - cdf + 1 >= 0.5)]
  67 serial_med
  68 ```
  69
  70 Below is a graph of the above results, containing the serial distribution as
  71 well as its mean and median, which could be used as estimates of the serial
  72 interval:
  73
  74 ```{r, dpi = 192, echo = FALSE}
  75 par(mar = c(4.1, 4.1, 0.5, 0.5))
  76
  77 # Serial distribution.
  78 plot(estimate$supp, estimate$pmf, xlab = "x", ylab = "p(x)",
  79   col = "black", lty = 1, type = "l"
  80 )
  81
  82 # Serial mean.
  83 abline(v = serial_mean, col = "blue", lty = 2)
  84 # Serial median.
  85 abline(v = serial_med, col = "red", lty = 2)
  86
  87 legend("topright",
  88   legend = c("Serial distribution", "Serial mean", "Serial median"),
  89   col = c("black", "blue", "red"),
  90   lty = c(1, 2, 2), cex = 0.5
  91 )
  92 ```