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

   1 ---
   2 title: "Sequential Bayes: Utilizing the Posterior Distribution"
   3 output: rmarkdown::html_vignette
   4 vignette: >
   5   %\VignetteIndexEntry{Sequential Bayes: Utilizing the Posterior 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 In the Sequential Bayes method, the probability distribution of R0 is updated
  22 sequentially from one case count to the next, starting from a (discretized)
  23 uniform prior. By default, the function `seq_bayes()` returns the mean of the
  24 last updated posterior distribution as its estimate of R0. However, by setting
  25 the parameter `post` to `TRUE`, it is possible to return the final distribution
  26 itself:
  27
  28 ```{r}
  29 # Daily case counts.
  30 cases <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)
  31
  32 posterior <- seq_bayes(cases, mu = 8, kappa = 7, post = TRUE)
  33 ```
  34
  35 First, the distribution can be used to retrieve the original estimate (had
  36 `post` been left to its default value of `FALSE`) by calculating its mean:
  37
  38 ```{r}
  39 # `supp` is the support of the distribution, and `pmf` is its probability mass
  40 # function.
  41 post_mean <- sum(posterior$supp * posterior$pmf)
  42 post_mean
  43
  44 # Verify that the following is true:
  45 post_mean == seq_bayes(cases, mu = 8, kappa = 7)
  46 ```
  47
  48 Another use of the posterior is to obtain an alternative estimate of R0. For
  49 instance, the following extracts the posterior mode rather than the mean:
  50
  51 ```{r}
  52 post_mode <- posterior$supp[which.max(posterior$pmf)]
  53 post_mode
  54 ```
  55
  56 Returning the posterior is suitable for visualization purposes. Below is a graph
  57 containing the uniform prior, final posterior distribution, posterior mean and
  58 posterior mode:
  59
  60 ```{r, dpi = 192, echo = FALSE}
  61 par(mar = c(4.1, 4.1, 0.5, 0.5))
  62
  63 # Posterior.
  64 plot(posterior$supp, posterior$pmf, xlab = "x", ylab = "p(x)",
  65   col = "black", lty = 1, type = "l"
  66 )
  67 # Uniform prior.
  68 segments(x0 = 0, x1 = 7, y0 = 1 / (7 / 0.01 + 1), y1 = 1 / (7 / 0.01 + 1),
  69   col = "orange"
  70 )
  71 # Posterior mean.
  72 abline(v = post_mean, col = "blue", lty = 2)
  73 # Posterior mode.
  74 abline(v = post_mode, col = "red", lty = 2)
  75
  76 legend("topright",
  77   legend = c("Prior", "Posterior", "Posterior mean", "Posterior mode"),
  78   col = c("orange", "black", "blue", "red"),
  79   lty = c(1, 1, 2, 2), cex = 0.5
  80 )
  81 ```