add sample_posterior_R function

NAMESPACE

@@ -13,6 +13,7 @@ export(init_mcmc_params)
 export(make_config)
 export(make_mcmc_control)
 export(overall_infectivity)
+export(sample_posterior_R)
 export(wallinga_teunis)
 import(grid)
 import(gridExtra)

NEWS.md

@@ -1,5 +1,10 @@
 # EpiEstim 2.2-0
 
+## NEW FUNCTIONS
+
+* `sample_posterior_R()` samples values of R from the posterior distribution of
+  an `estimate_R` object (#70, @acori)
+
 ## MISC
 
 * Added a `NEWS.md` file to track changes to the package. (#74, @zkamvar)

R/sample_posterior_R.R

@@ -0,0 +1,56 @@
+#' sample from the posterior R distribution
+#'
+#' @param R an \code{estimate_R} object from the estimate_r function
+#'   function.
+#' 
+#' @param n an integer specifying the number of samples to be taken from the
+#'   gamma distribution.
+#' 
+#' @param window an integer (or sequence of integers) specifying the window(s) 
+#'   from which to estimate R. Defaults to the first window. If multiple windows
+#'   are specified, the resulting samples will be drawn from several 
+#'   distributions.
+#' 
+#' @return n values of R from the posterior R distribution
+#'
+#' @author Anne Cori
+#' @export
+#'
+#' @examples
+#'
+#'
+#' ## load data on pandemic flu in a school in 2009
+#' data("Flu2009")
+#'
+#' ## estimate the reproduction number (method "non_parametric_si")
+#' ## when not specifying t_start and t_end in config, they are set to estimate
+#' ## the reproduction number on sliding weekly windows                          
+#' res <- estimate_R(incid = Flu2009$incidence, 
+#'                   method = "non_parametric_si",
+#'                   config = make_config(list(si_distr = Flu2009$si_distr)))
+#' 
+#' ## Sample R from the first weekly window
+#' win <- 1L
+#' R_median <- res$R$`Median(R)`[win]
+#' R_CrI <- c(res$R$`Quantile.0.025(R)`[win], res$R$`Quantile.0.975(R)`[win])
+#'
+#' set.seed(2019-06-06) # fixing the random seed for reproducibility
+#' R_sample <- sample_posterior_R(res, n = 1000, window = win)
+#' hist(R_sample, col = "grey", main = "R sampled from the first weekly window")
+#' abline(v = R_median, col = "red")       # show the median estimated R
+#' abline(v = R_CrI, col = "red", lty = 2) # show the 95%CrI of R
+sample_posterior_R <- function(R, n = 1000, window = 1L) {
+
+  if (!inherits(R, "estimate_R")) {
+    stop("R must be generated from the estimate_r function.")
+  }
+  mu    <- R$R$`Mean(R)`[window]
+  sigma <- R$R$`Std(R)`[window]
+  # Gamma distribution needs a shape and scale, which are transformations
+  # of the mean (mu) and standard deviation (sigma) 
+  cv    <- sigma / mu # coefficient of variation
+  shape <- 1 / (cv ^ 2)
+  scale <- mu * cv ^ 2
+  rgamma(n, shape = shape, scale = scale)
+
+}

tests/testthat/test-stample_posterior_R.R

@@ -0,0 +1,44 @@
+
+data("Flu2009")
+
+## estimate the reproduction number (method "non_parametric_si")
+## when not specifying t_start and t_end in config, they are set to estimate
+## the reproduction number on sliding weekly windows                          
+res <- estimate_R(incid = Flu2009$incidence, 
+                  method = "non_parametric_si",
+                  config = make_config(list(si_distr = Flu2009$si_distr)))
+
+
+test_that("only estimate_R objects are accepted", {
+
+  expect_error(sample_posterior_R(res$R), "R must be generated from the estimate_r function")                  
+
+})
+
+test_that("sampling returns n samples of R", {
+
+  set.seed(2019)            
+  R_sample <- sample_posterior_R(res, n = 1001, window = 1L)
+  expect_length(R_sample, 1001)
+
+})
+
+
+test_that("different windows can be specified", {
+
+  set.seed(2019)
+  R_sample1 <- sample_posterior_R(res, n = 1001, window = 1L)
+
+  set.seed(2019)
+  R_sample2 <- sample_posterior_R(res, n = 1001, window = 1L)
+
+  set.seed(2019)
+  R_sample3 <- sample_posterior_R(res, n = 1001, window = 21L)
+
+  # Samples from the same window should be the same
+  expect_equal(R_sample1, R_sample2)
+
+  # Samples from a different window should be different
+  expect_failure(expect_equal(R_sample1, R_sample3))
+
+})