From 59de2553220ffff0a62ed3e454876df2bb93916b Mon Sep 17 00:00:00 2001 From: Naeem Model Date: Sat, 3 Feb 2024 15:17:20 +0000 Subject: Move seqB posterior example to vignettes --- R/seq_bayes.R | 3 ++ man/seq_bayes.Rd | 4 +++ vignettes/seq_bayes_post.Rmd | 80 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 87 insertions(+) create mode 100644 vignettes/seq_bayes_post.Rmd diff --git a/R/seq_bayes.R b/R/seq_bayes.R index 107d428..2f87121 100644 --- a/R/seq_bayes.R +++ b/R/seq_bayes.R @@ -55,6 +55,9 @@ #' [Bettencourt and Riberio (PloS One, 2008)]( #' https://doi.org/10.1371/journal.pone.0002185) #' +#' @seealso `vignette("seq_bayes_post", package = "Rnaught")` for examples of +#' using the posterior distribution. +#' #' @export #' #' @examples diff --git a/man/seq_bayes.Rd b/man/seq_bayes.Rd index 41445a9..8dc824c 100644 --- a/man/seq_bayes.Rd +++ b/man/seq_bayes.Rd @@ -90,3 +90,7 @@ estimate == sum(posterior$supp * posterior$pmf) \references{ \href{https://doi.org/10.1371/journal.pone.0002185}{Bettencourt and Riberio (PloS One, 2008)} } +\seealso{ +\code{vignette("seq_bayes_post", package = "Rnaught")} for examples of +using the posterior distribution. +} diff --git a/vignettes/seq_bayes_post.Rmd b/vignettes/seq_bayes_post.Rmd new file mode 100644 index 0000000..7c614ba --- /dev/null +++ b/vignettes/seq_bayes_post.Rmd @@ -0,0 +1,80 @@ +--- +title: "Sequential Bayes: Utilizing the Posterior Distribution" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Sequential Bayes: Utilizing the Posterior Distribution} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup, include = FALSE} +library(Rnaught) +``` + +In the Sequential Bayes method, the probability distribution of R0 is updated +sequentially from one case count to the next, starting from a (discretized) uniform prior. By +default, the function `seq_bayes` returns the mean of the last updated posterior +distribution as its estimate of R0. However, by setting the parameter `post` to +`TRUE`, it is possible to return the final distribution itself: + +```{r} +# Daily case counts. +cases <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) + +posterior <- seq_bayes(cases, mu = 8, kappa = 7, post = TRUE) +``` + +First, the distribution can be used to retrieve the original estimate (had +`post` been left to its default value of `FALSE`) by calculating its mean: + +```{r} +# `supp` is the support of the distribution, and `pmf` is its probability mass +# function. +post_mean <- sum(posterior$supp * posterior$pmf) +post_mean + +# Verify that the following is true: +post_mean == seq_bayes(cases, mu = 8, kappa = 7) +``` + +Another use of the posterior is to obtain an alternative estimate of R0. For +instance, the following extracts the posterior mode rather than the mean: + +```{r} +post_mode <- posterior$supp[which.max(posterior$pmf)] +post_mode +``` + +Returning the posterior is suitable for visualization purposes. Below is a graph +containing the uniform prior, final posterior distribution, posterior mean and +posterior mode: + +```{r, dpi = 192, echo = FALSE} +par(mar = c(4.1, 4.1, 0.5, 0.5)) + +# Posterior. +plot(posterior$supp, posterior$pmf, xlab = "x", ylab = "p(x)", + col = "black", lty = 1, type = "l" +) +# Uniform prior. +segments(x0 = 0, x1 = 7, y0 = 1 / (7 / 0.01 + 1), y1 = 1 / (7 / 0.01 + 1), + col = "orange" +) +# Posterior mean. +abline(v = post_mean, col = "blue", lty = 2) +# Posterior mode. +abline(v = post_mode, col = "red", lty = 2) + +legend("topright", + legend = c("Prior", "Posterior", "Posterior mean", "Posterior mode"), + col = c("orange", "black", "blue", "red"), + lty = c(1, 1, 2, 2), cex = 0.5 +) +``` -- cgit v1.2.3