fix EpiNow2::forecast_secondary() and EpiNow2::plot.estimate_secondary()

episodes/create-forecast.Rmd

@@ -298,22 +298,23 @@ The final key input is the delay distribution between the primary and secondary
 There are further function inputs to `estimate_secondary()` which can be specified, including adding an observation process, see `?EpiNow2::estimate_secondary` for detail on these options. 
 
 To find the model fit between cases and deaths: 
-```{r}
+```{r,warning=FALSE,message=FALSE}
 # Estimate from day 31 to day 60 of this data
 cases_to_estimate <- reported_cases_deaths %>%
  slice(31:60)
 
+# Delay distribution between case report and deaths
+delay_report_to_death <- EpiNow2::Gamma(
+ mean = EpiNow2::Normal(mean = 14, sd = 0.5),
+ sd = EpiNow2::Normal(mean = 5, sd = 0.5),
+ max = 30
+)
+
 # Estimate secondary cases
 estimate_cases_to_deaths <- EpiNow2::estimate_secondary(
  data = cases_to_estimate,
  secondary = EpiNow2::secondary_opts(type = "incidence"),
- delays = EpiNow2::delay_opts(
- EpiNow2::Gamma(
- mean = 14,
- sd = 5,
- max = 30
- )
- )
+ delays = EpiNow2::delay_opts(delay_report_to_death)
 )
 ```
 
@@ -339,12 +340,12 @@ To use this model fit to forecast deaths, we pass a data frame consisting of the
 # Forecast from day 61 to day 90
 cases_to_forecast <- reported_cases_deaths %>%
  dplyr::slice(61:90) %>%
- dplyr::select(date, value = primary)
+ dplyr::mutate(value = primary)
 ```
 
 To forecast, we use the model fit `estimate_cases_to_deaths`: 
 
-```{r}
+```{r,warning=FALSE,message=FALSE}
 # Forecast secondary cases
 deaths_forecast <- EpiNow2::forecast_secondary(
  estimate = estimate_cases_to_deaths,
@@ -354,6 +355,34 @@ deaths_forecast <- EpiNow2::forecast_secondary(
 plot(deaths_forecast)
 ```
 
+:::::::::::::::: spoiler
+
+### make a forecast plot
+
+```{r}
+deaths_forecast %>%
+ purrr::pluck("predictions") %>%
+ ggplot(aes(x = date, y = secondary)) +
+ geom_col(
+ fill = "grey", col = "white",
+ show.legend = FALSE, na.rm = TRUE
+ ) +
+ geom_ribbon(aes(ymin = lower_90, ymax = upper_90),
+ alpha = 0.2, linewidth = 1) +
+ geom_ribbon(aes(ymin = lower_50, ymax = upper_50),
+ alpha = 0.4, linewidth = 1) +
+ geom_ribbon(aes(ymin = lower_20, ymax = upper_20),
+ alpha = 0.6, linewidth = 1) +
+ theme_bw() +
+ labs(y = "Reports per day", x = "Date") +
+ scale_x_date(date_breaks = "week", date_labels = "%b %d") +
+ scale_y_continuous(labels = scales::comma) +
+ theme(axis.text.x = ggplot2::element_text(angle = 90))
+```
+
+
+::::::::::::::::
+
 The plot shows the forecast secondary observations (deaths) over the dates which we have recorded cases for. 
 It is also possible to forecast deaths using forecast cases, here you would specify `primary` as the `estimates` output from `estimate_infections()`.