Create vignette for WP serial distribution

2 files changed, 95 insertions, 0 deletions

diff --git a/R/wp.R b/R/wp.R
index 674a723..3448911 100644
--- a/R/wp.R
+++ b/R/wp.R
@@ -75,6 +75,9 @@
 #' [White and Pagano (Statistics in Medicine, 2008)](
 #' https://doi.org/10.1002/sim.3136)
 #'
+#' @seealso `vignette("wp_serial", package="Rnaught")` for examples of using the
+#'   serial distribution.
+#'
 #' @importFrom stats pgamma qgamma
 #'
 #' @export
diff --git a/vignettes/wp_serial.Rmd b/vignettes/wp_serial.Rmd
new file mode 100644
index 0000000..63f89a7
--- /dev/null
+++ b/vignettes/wp_serial.Rmd
@@ -0,0 +1,92 @@
+---
+title: "White and Pagano: Utilizing the Serial Distribution"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{White and Pagano: Utilizing the Serial Distribution}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup, include = FALSE}
+library(Rnaught)
+```
+
+The serial distribution of an infectious disease is the distribution of the time
+from when an infectious individual -- the infector -- becomes symptomatic, to
+when another individual who is infected by the infector becomes symptomatic. The
+serial interval refers to a range of likely values from this distribution,
+although it is typically reported as the mean.
+
+In the White and Pagano method, the serial distribution is assumed to be a
+discretized, finite version of a gamma distribution. Setting the parameter
+`serial` to `TRUE` causes this discretized distribution to be returned in
+addition to the estimate of R0. Furthermore, the method can be used whether or
+not the serial interval (specified as the parameter `mu`) is known. When `mu` is
+specified, it is taken to be the mean of a continuous gamma distribution (i.e.,
+before the discretization). As such, the mean computed from the returned serial
+distribution may differ slightly from `mu`:
+
+```{r}
+# Case counts.
+cases <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
+
+estimate <- wp(cases, mu = 3.333, serial = TRUE)
+
+# `supp` is the support of the distribution, and `pmf` is its probability mass
+# function.
+sum(estimate$supp * estimate$pmf)
+```
+
+When `mu` is unspecified (left to its default value of `NA`), the method
+performs a maximum likelihood estimation over all (discretized) gamma
+distributions via a grid search, whose range of parameters are specified via
+`grid_length`, `max_shape` and `max_scale` (see `?wp` for more details). It is
+useful to return the estimated serial distribution in this case, as it can
+provide estimates of the serial interval when it is unknown:
+
+```{r}
+# The grid search parameters specified below are the default values.
+estimate <- wp(cases, serial = TRUE,
+  grid_length = 100, max_shape = 10, max_scale = 10
+)
+
+serial_mean <- sum(estimate$supp * estimate$pmf)
+serial_mean
+
+# Compute the (discrete) median for an alternative estimate of the serial
+# interval.
+cdf <- cumsum(estimate$pmf)
+serial_med <- estimate$supp[which(cdf >= 0.5 & estimate$pmf - cdf + 1 >= 0.5)]
+serial_med
+```
+
+Below is a graph of the above results, containing the serial distribution as
+well as its mean and median, which could be used as estimates of the serial
+interval:
+
+```{r, dpi = 192, echo = FALSE}
+par(mar = c(4.1, 4.1, 0.5, 0.5))
+
+# Serial distribution.
+plot(estimate$supp, estimate$pmf, xlab = "x", ylab = "p(x)",
+  col = "black", lty = 1, type = "l"
+)
+
+# Serial mean.
+abline(v = serial_mean, col = "blue", lty = 2)
+# Serial median.
+abline(v = serial_med, col = "red", lty = 2)
+
+legend("topright",
+  legend = c("Serial distribution", "Serial mean", "Serial median"),
+  col = c("black", "blue", "red"),
+  lty = c(1, 2, 2), cex = 0.5
+)
+```