Updating transmissibility template

inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/EpiEstim.Rmd

@@ -10,6 +10,15 @@ library(EpiEstim)
 
 ### Explanations
 
+The package *EpiEstim* is used to estimate the effective, time-varying reproduction number, $Rt$,
+using a branching process method that requires daily incidence data and an estimate of the disease's serial interval. Through this method, $Rt$ is estimated over weekly sliding windows, i.e., successive overlapping 7-day periods. $Rt$ estimates can be biased by reporting delays. For this reason, this report excludes the last `r params$incomplete_days` days included in the incidence data as incomplete, and retains a total of `r params$r_estim_window` days to obtain the most recent estimate of the time varying reproduction number.
+
+
+$Rt$ is defined as the average number of secondary cases that an infected individual would infect if
+conditions remained as they were at time _t_. It represents the transmission potential of a disease at a specified timepoint. Observing changes in $Rt$ provides valuable information about variations in transmissibility during an epidemic, and can be used to assess the effectiveness of control interventions aimed at reducing virus transmission. 
+For more information on *Rt* and its interpretation, see [this paper by the Royal Society](https://royalsociety.org/-/media/policy/projects/set-c/set-covid-19-R-estimates.pdf).
+
+
 ### Results
 
 ```{r rt-wrappers-prep}
@@ -69,6 +78,8 @@ associated 95% credibility intervals
 * a table summarising these results for the last `r params$r_estim_window` days
 
 ```{r rt-estim-global}
+# This code uses the {EpiEstim} function `estimate_R` to generate an estimate of
+# the global Rt over time
 res_epiestim_global <- dat_i_day %>%
   regroup() %>%
   get_count_value() %>%
@@ -77,7 +88,7 @@ res_epiestim_global <- dat_i_day %>%
 ```
 
 ```{r rt-plot-global}
-# Graph of all values over time
+# Plot to visualise the global Rt over time
 ggplot(res_epiestim_global, aes(x = end)) +
   geom_ribbon(aes(ymin = lower, ymax = upper), fill = pale_green) +
   geom_line(aes(y = median), color = dark_green) +
@@ -91,7 +102,7 @@ ggplot(res_epiestim_global, aes(x = end)) +
 ```
 
 ```{r rt-table-global}
-# Table
+# This code generates a table with Rt values on the last params$r_estim_window
 res_epiestim_global %>%
   tail(params$r_estim_window) %>%
   mutate(
@@ -116,22 +127,31 @@ res_epiestim_global %>%
 
 #### Transmissibility by group
 
+These analyses present Rt results for the stratified incidence. Results include:
+
+* a graph of the estimates of $R_t$ of the observed data over time per group, with
+associated 95% credibility intervals
+
+* a plot with the latest $R_t$ estimate per group, with associated 95% credibility intervals
+
+* a table with the latest $R_t$ estimates per group, with associated 95% credibility intervals
+
 ```{r rt-estim-group}
-# Get results by group, keeping only the last 14 days of data
+# This code estimates Rt by group_var, keeping only the last 14 days of data
 res_epiestim_group <- dat_i_day %>%
   nest(data = c(get_date_index_name(.), get_count_value_name(.))) %>%
   mutate(
     res_epiestim = map(data, ~ wrap_res(
       estimate_R(.x$count, config = ee_config),
-      dat_i_day)
-    )
+      dat_i_day
+    ))
   ) %>%
   unnest(res_epiestim) %>%
   dplyr::select(-data)
 ```
 
 ```{r rt-plot-group, fig.height = 5 / 3 * n_groups}
-# Graph of all values over time
+# Plot of Rt estimates for each group_var over time
 ggplot(res_epiestim_group, aes(x = end)) +
   geom_ribbon(aes(ymin = lower, ymax = upper), fill = pale_green) +
   geom_line(aes(y = median), color = dark_green) +
@@ -146,25 +166,26 @@ ggplot(res_epiestim_group, aes(x = end)) +
 ```
 
 ```{r rt-plot-estimates-group}
-# Plot of latest estimates
+# Plot of latest estimates of Rt for each group_var
 res_epiestim_group %>%
   filter(end == max(end)) %>%
   ggplot(aes(y = .data[[group_var]]), fill = custom_grey) +
-    geom_point(aes(x = median), color = dark_green) +
-    geom_errorbar(aes(xmin = lower, xmax = upper), color = dark_green) +
-    geom_vline(xintercept = 1, color = dark_pink) +
-    labs(
-      title = "Latest estimates of Rt",
-      subtitle = sprintf(
-        "As of %s",
-        format(max(get_dates(dat_i_day)), "%d %B %Y")
-      ),
-      y = "",
-      x = "Instantaneous Reproduction Number (Rt)"
-    )
+  geom_point(aes(x = median), color = dark_green) +
+  geom_errorbar(aes(xmin = lower, xmax = upper), color = dark_green) +
+  geom_vline(xintercept = 1, color = dark_pink) +
+  labs(
+    title = "Latest estimates of Rt",
+    subtitle = sprintf(
+      "As of %s",
+      format(max(get_dates(dat_i_day)), "%d %B %Y")
+    ),
+    y = "",
+    x = "Instantaneous Reproduction Number (Rt)"
+  )
 ```
 
 ```{r rt-table-group}
+# This code generates a table with the latest Rt estimates for each group_var
 res_epiestim_group %>%
   filter(end == max(end)) %>%
   mutate(
@@ -176,7 +197,7 @@ res_epiestim_group %>%
       upper
     )
   ) %>%
-  dplyr::select(-c(sd, lower, upper)) %>%
+  dplyr::select(-c(sd, lower, upper, count_variable)) %>%
   rename(
     "mean $R$" = mean,
     "median $R$" = median

inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/EpiNow2.Rmd

@@ -10,28 +10,20 @@ library(EpiNow2)
 
 ### Explanations
 
+The package *EpiNow2* is used to estimate the effective, time-varying reproduction number, $Rt$, using a Bayesian inference method that requires data on observed infections, and an estimate of the generation time, which is approximated using the disease's serial interval distribution. $Rt$ estimates can be biased by reporting delays. For this reason, this report excludes the last `r params$incomplete_days` days included in the incidence data as incomplete,
+and retains a total of `r params$r_estim_window` days to obtain the most recent estimate of the time varying reproduction number.
+
+$Rt$ is defined as the average number of secondary cases that an infected individual would infect if
+conditions remained as they were at time _t_. It represents the transmission potential of a disease at a specified timepoint. Observing changes in $Rt$ provides valuable information about variations in transmissibility during an epidemic, and can be used to assess the effectiveness of control interventions aimed at reducing virus transmission. 
+For more information on *Rt* and its interpretation, see [this paper by the Royal Society](https://royalsociety.org/-/media/policy/projects/set-c/set-covid-19-R-estimates.pdf).
+
 ### Results
 
 ```{r rt-wrappers-prep}
-# Approximate serial interval with gamma distribution since this is what EpiNow2
-# will use for generation_time
-si_gamma <- epiparameter::extract_param(
-  type = "range",
-  values = c(median(si$prob_dist$r(1e3)),
-             min(si$prob_dist$r(1e3)),
-             max(si$prob_dist$r(1e3))),
-  distribution = params$si_dist,
-  samples = 1e3
-)
-si_gamma <- epiparameter::convert_params_to_summary_stats(
-  distribution = "gamma", shape = si_gamma[[1]], scale = si_gamma[[2]])
-
-generation_time <- list(
-  mean = si_gamma$mean,
-  mean_sd = 0,
-  sd = si_gamma$sd,
-  sd_sd = 0,
-  max = max(si_x)
+# This code takes the serial interval distribution of params$disease_name and
+# uses it to approximate the disease's generation time
+generation_time <-  generation_time_opts(
+  dist_spec(pmf = density(si, at = si_x))
 )
 ```
 
@@ -43,16 +35,19 @@ stratification. Results include:
 * a graph of the estimates of $R_t$ of the observed data over time, with
 associated 95% credibility intervals
 
-* a table summarising these results for the last `r params$r_estim_window` days
+* a table summarising $Rt$ values for the last `r params$r_estim_window` days of available data
 
 ```{r rt-estim-global}
+# This code uses {EpiNow2}'s function `epinow()` to estimate Rt over time
+# from daily incidence data
 # Warning, this step will take 2 cores of your computer and may be running
 # for a long time!
 res_epinow2_global <- dat_i_day %>%
   regroup() %>%
-  dplyr::rename(
-    confirm = .data[[count_var]],
-    date    = date_index
+  dplyr::mutate(
+    confirm = count,
+    date    = as.Date(date_index),
+    .keep   = "unused"
   ) %>%
   epinow(
     generation_time = generation_time,
@@ -67,6 +62,7 @@ res_epinow2_global <- dat_i_day %>%
 ```
 
 ```{r rt-plot-global}
+# Plot of Rt values over time for global incidence data (without stratification)
 plot(res_epinow2_global, "R") +
   scale_x_date(breaks = scales::pretty_breaks(n = 8), date_label = "%d %b %Y") +
   theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
@@ -78,6 +74,8 @@ plot(res_epinow2_global, "R") +
 ```
 
 ```{r rt-table-global}
+# This code generates a table that contains Rt estimates for the last
+# params$r_estim_window
 res_epinow2_global$estimates$summarised %>%
   dplyr::filter(variable == "R") %>%
   tail(params$r_estim_window) %>%
@@ -98,15 +96,27 @@ res_epinow2_global$estimates$summarised %>%
 
 #### Transmissibility by group
 
+These analyses present results for stratified incidence, i.e., by `r group_var`. Results include:
+
+* a graph of the estimates of $R_t$ of the observed data over time, with associated 95% credibility intervals, stratified by `r group_var`
+
+* a table with $Rt$ estimates for each `r group_var`, excluding the most recent `r params$incomplete_days` days from the available data 
+
 ```{r rt-estim-group, message = FALSE}
+# This code uses {EpiNow2}'s function `epinow()` to estimate Rt over time
+# from daily incidence data, stratified by group_var
+# Warning, this step will take 2 cores of your computer and may be running
+# for a long time!
 res_epinow2_group <- dat_i_day %>%
+  select("date_index", "region", "count") %>%
   dplyr::rename(
-    confirm = .data[[count_var]],
+    confirm = .data[["count"]],
     date    = date_index,
     region  = .data[[group_var]]
   ) %>%
+  dplyr::mutate(date = as.Date(date)) %>%
   regional_epinow(
-    generation_time = generation_time,
+    generation_time = generation_time_opts(generation_time),
     stan = stan_opts(samples = 1e3, chains = 2),
     rt = rt_opts(gp_on = "R0"),
     horizon = 0,
@@ -118,14 +128,12 @@ res_epinow2_group <- dat_i_day %>%
 ```
 
 ```{r rt-plot-group}
+# Plot of Rt estimates over time, stratified by group_var
 res_epinow2_group$summary$plots$R
 ```
 
-```{r rt-plot-estimates-group}
-# Missing for now
-```
-
 ```{r rt-table-group}
+# This code generates a table with an overall estimate of Rt per group_var
 res_epinow2_group$summary$summarised_results$data %>%
   dplyr::filter(metric == "Effective reproduction no.") %>%
   dplyr::select(region, mean, median, lower_95, upper_95) %>%

inst/rmarkdown/templates/transmissibility/skeleton/rmdchunks/R0.Rmd

@@ -10,9 +10,17 @@ library(R0)
 
 ### Explanations
 
+The package *R0* is used to estimate the effective, time-varying reproduction number, $Rt$, implementing a time-dependent method, described by [Wallinga and Teunis](https://academic.oup.com/aje/article/160/6/509/79472?login=false). This method requires incidence data and an estimate of the generation time, which here is approximated using the disease's serial interval distribution.
+$Rt$ estimates can be biased by reporting delays. For this reason, this report excludes the last `r params$incomplete_days` days included in the incidence data as incomplete.
+
+$Rt$ is defined as the average number of secondary cases that an infected individual would infect if
+conditions remained as they were at time _t_. It represents the transmission potential of a disease at a specified timepoint. Observing changes in $Rt$ provides valuable information about variations in transmissibility during an epidemic, and can be used to assess the effectiveness of control interventions aimed at reducing virus transmission. 
+For more information on *Rt* and its interpretation, see [this paper by the Royal Society](https://royalsociety.org/-/media/policy/projects/set-c/set-covid-19-R-estimates.pdf).
+
 ### Results
 
 ```{r rt-wrappers-prep}
+# This code uses the serial interval to approximate the generation time
 mGT <- generation.time("empirical", si$prob_dist$d(si_x))
 
 # Helper function to get R0 outputs in a tidy format, compatible with pipelines
@@ -41,19 +49,23 @@ associated 95% credibility intervals
 * a table summarising these results for the last `r params$r_estim_window` days
 
 ```{r rt-estim-global}
+# This code uses the {R0} function `estimate.R()` to generate Rt estimates over
+# time using daily incidence data (dat_i_day)
 res_r0_global <- dat_i_day %>%
   regroup() %>%
   mutate(r0_quantiles(estimate.R(
-      get_count_value(.),
-      mGT,
-      date_index,
-      begin = 1L, end = length(date_index),
-      nsim = 1e4, method = "TD"
-    )$estimates$TD)
+    get_count_value(.),
+    mGT,
+    date_index,
+    begin = 1L, end = length(date_index),
+    nsim = 1e4, method = "TD"
+  )$estimates$TD)
   )
 ```
 
 ```{r rt-plot-global}
+# Plot of global Rt estimates over time for daily incidence data included in
+# dat_i_day
 res_r0_global %>%
   na.omit() %>%
   ggplot(aes(x = date_index, y = median, ymin = lower, ymax = upper)) +
@@ -69,6 +81,8 @@ res_r0_global %>%
 ```
 
 ```{r rt-table-global}
+# This code generates a table with Rt estimates for the last
+# `r params$r_estim_window` days included in dat_i_day
 res_r0_global %>%
   tail(params$r_estim_window) %>%
   mutate(
@@ -85,13 +99,24 @@ res_r0_global %>%
     "median $R$" = median
   ) %>%
   set_names(toupper) %>%
-  kbl() %>%
+  kbl(row.names = FALSE) %>%
   kable_paper("striped", font_size = 18, full_width = FALSE)
 ```
 
 #### Transmissibility by group
 
+These analyses present results for incidence stratified by `r group_var`. Results in this section include:
+
+* a graph of the estimates of $R_t$ of the observed data over time for each `r group_var`, with
+associated 95% credibility intervals
+
+* a plot with the latest $R_t$ results by `r group_var`
+
+* a table with the latest $R_t$ estimates by `r group_var`
+
 ```{r rt-estim-group}
+# This code uses the {R0} function `estimate_R()` to obtain Rt estimates by
+# group_var over time
 res_r0_group <- dat_i_day %>%
   regroup(groups = get_group_names(.)) %>%
   arrange(date_index) %>%
@@ -108,6 +133,7 @@ res_r0_group <- dat_i_day %>%
 ```
 
 ```{r rt-plot-group}
+# Plot of Rt estimates over time by group_var
 res_r0_group %>%
   na.omit() %>%
   ggplot(aes(x = date_index, y = median, ymin = lower, ymax = upper)) +
@@ -124,26 +150,27 @@ res_r0_group %>%
 ```
 
 ```{r rt-plot-estimates-group}
-# Plot of latest estimates
+# Plot of latest estimates of Rt by group_var
 res_r0_group %>%
   na.omit() %>%
   filter(date_index == max(date_index)) %>%
   ggplot(aes(y = .data[[group_var]]), fill = custom_grey) +
-    geom_point(aes(x = median), color = dark_green) +
-    geom_errorbar(aes(xmin = lower, xmax = upper), color = dark_green) +
-    geom_vline(xintercept = 1, color = dark_pink) +
-    labs(
-      title = "Latest estimates of Rt",
-      subtitle = sprintf(
-        "As of %s",
-        format(max(get_dates(dat_i_day)), "%d %B %Y")
-      ),
-      y = "",
-      x = "Instantaneous Reproduction Number (Rt)"
-    )
+  geom_point(aes(x = median), color = dark_green) +
+  geom_errorbar(aes(xmin = lower, xmax = upper), color = dark_green) +
+  geom_vline(xintercept = 1, color = dark_pink) +
+  labs(
+    title = "Latest estimates of Rt",
+    subtitle = sprintf(
+      "As of %s",
+      format(max(get_dates(dat_i_day)), "%d %B %Y")
+    ),
+    y = "",
+    x = "Instantaneous Reproduction Number (Rt)"
+  )
 ```
 
 ```{r rt-table-group}
+# This code generates a table with latest estimates of Rt by group_var
 res_r0_group %>%
   na.omit() %>%
   filter(date_index == max(date_index)) %>%