EpiNow2 icon indicating copy to clipboard operation
EpiNow2 copied to clipboard

error in `plot.estimate_secondary()` to plot the output from `forecast_secondary()`

Open avallecam opened this issue 8 months ago • 0 comments

Summary: Error in plot.estimate_secondary() to plot the output from forecast_secondary() due to filtering or subsetting action in the first line of the plotting function.

Description:

This undesired behaviour seems to be visible in the estimate_secondary() reference manual example line: plot(prev_preds, new_obs = cases, from = "2020-06-01")

In the most recent version, the [!is.na(secondary)] was added, which drops all the forecast rows in the dataset.

https://github.com/epiforecasts/EpiNow2/blob/0b776fd1457c13cafd936d760e1eb045570cbd3e/R/estimate_secondary.R#L350-L355

a previous commit did not have this additional subsetting.

https://github.com/epiforecasts/EpiNow2/blob/df1fdc8bcb361da308037fbb9c8deb6f3fc45235/R/estimate_secondary.R#L393-L395

Reproducible Steps:


library(EpiNow2)
#> 
#> Attaching package: 'EpiNow2'
#> The following object is masked from 'package:stats':
#> 
#>     Gamma
library(cfr)
library(tidyverse)

# data --------------------------------------------------------------------

reported_cases_deaths <- cfr::ebola1976 %>% 
  dplyr::rename(primary = cases, secondary = deaths)

# Estimate from day 1 to day 35 of this data
cases_to_estimate <- reported_cases_deaths %>%
  slice(1:35)

# Forecast from day 36 to day 73
cases_to_forecast <- reported_cases_deaths %>%
  dplyr::slice(36:73) %>% 
  dplyr::mutate(value = primary)


# EpiNow2 -----------------------------------------------------------------

# delay
delay_uncertain <- EpiNow2::Gamma(
  mean = EpiNow2::Normal(mean = 14, sd = 0.5),
  sd = EpiNow2::Normal(mean = 5, sd = 0.5),
  max = 30
)
#> Warning in new_dist_spec(params, "gamma"): Uncertain gamma distribution
#> specified in terms of parameters that are not the "natural" parameters of the
#> distribution (shape, rate). Converting using a crude and very approximate
#> method that is likely to produce biased results. If possible, it is preferable
#> to specify the distribution directly in terms of the natural parameters.

# estimate secondary cases
secondary_uncertain <- EpiNow2::estimate_secondary(
  data = cases_to_estimate,
  delays = EpiNow2::delay_opts(delay_uncertain)
)
#> WARN [2024-06-22 02:05:03] estimate_secondary (chain: 1): Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess -

# forecast secondary cases
forecast_uncertain <- EpiNow2::forecast_secondary(
  estimate = secondary_uncertain,
  primary = cases_to_forecast
)


# ISSUE: plot() does not plot the forecast --------------------------------

# only plot the observed cases
plot(forecast_uncertain)
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf



# possible solution? ------------------------------------------------------

# from
# ?EpiNow2:::plot.estimate_secondary()

x = forecast_uncertain
predictions <- data.table::copy(x$predictions)#[!is.na(secondary)]

library(ggplot2)

plot <- predictions %>% 
  ggplot(aes(x = date, y = secondary)) + 
  geom_col(
    fill = "grey", col = "white",
    show.legend = FALSE, na.rm = TRUE
  )

plot <- EpiNow2:::plot_CrIs(plot, extract_CrIs(predictions), alpha = 0.6, 
                  linewidth = 1)

plot +
  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))

Created on 2024-06-22 with reprex v2.1.0

EpiNow2 Version:

packageVersion(pkg = "EpiNow2")
#> [1] '1.5.2'

avallecam avatar Jun 22 '24 01:06 avallecam