non-continuity between the estimated and forecast growth rate credible intervals

[avallecam]

Summary:
CRAN version 1.4.0 and CRAN version 1.5.2 generate different figures from growth rate estimates.

Description:

• LEFT: from previous CRAN version 1.4.0 (generated in this commit md-output)
• RIGHT: current CRAN version 1.5.2 (generated in this commit md-output)

• LEFT: the credible intervale in the orange "estimate based on partial data" range increases towards Apr 27, and then the Forecast extends this uncertainty in the same magnitude range.
• RIGHT: the credible intervals in the same orange range have a higher range, but the Forecast does not seem to follow the uncertainty of the previous days.

This is currently visible in the tutorial episode but not in the reference manual.

Reproducible Steps:

With the code below, you can reproduce the figure.

library(EpiNow2)
library(cfr)
library(tidyverse)

withr::local_options(list(mc.cores = 4))

## generation time ------------------------------------------------

covid_serialinterval_epinow <- EpiNow2::LogNormal(
  meanlog = 1.39,
  sdlog = 0.57,
  max = 15
)

covid_serialinterval_epinow

# data 01 - 1:60 -----------------------------------------------------------------

data01_time02 <- EpiNow2::example_confirmed %>% 
  dplyr::slice(1:60)

## 00 delays -------------------------------------------------------

data01_time02_delays00 <- EpiNow2::epinow(
  # cases
  data = data01_time02,
  # delays
  generation_time = EpiNow2::generation_time_opts(covid_serialinterval_epinow),
  # stan 
  stan = EpiNow2::stan_opts(samples = 1000, chains = 2)
)

plot(data01_time02_delays00$plots$growth_rate) 

EpiNow2 Version:

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

---

[seabbs]

I can see a bug here but first are you looking at fits with the same seed in both R and stan?

[avallecam]

by looking at the code, possibly not, but these look similar when comparing different renderings of the same version. For R seed, I use set.seed(); how do you suggest setting the seed in stan?

[sbfnk]

Am I understanding correctly that you're getting different outputs when rendering repeatedly in spite of using set.seed()? That shouldn't be the case (and if it is then it's a bug) - the R seed is passed on to stan.

[avallecam]

Am I understanding correctly that you're getting different outputs when rendering repeatedly in spite of using set.seed()?

No, this is not the case, I just checked it using the reprex above and the output is stable without specifying any seed.

[sbfnk]

I think this might be an effect of the changed growth rate calculation in #610

The new version is the better version I think but it does mean that the non-smooth change in R translates into a discontinuity in growth (not 100% sure though - might be worth testing in simulation).

@avallecam can you check what it looks like when using rt = rt_opts(future = "project")?

[sbfnk]

might be worth testing in simulation

This is what we're seeing I think:

library("ggplot2")
library("dplyr")
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library("tidyr")

source("https://raw.githubusercontent.com/nfidd/nfidd/main/functions/renewal.r")
source("https://raw.githubusercontent.com/nfidd/nfidd/main/functions/censored-delay-pmf.r")


gen_pmf <- censored_delay_pmf(rgen = rgamma, max = 14, shape = 5, rate = 1)
R <- c(rep(1, 20), 1 + seq(0, 1, by = 0.05)**2, rep(2, 20))
I <- renewal(I0 = 1, R = R, gen_time = gen_pmf)
logI <- log(I)
## discard initial bit
growth <- c(NA_real_, diff(logI))

df <- tibble(t = seq_along(I), logI = logI, I = I, R = R, growth = growth) |>
  filter(t > 20) |> ## discard initial bit
  pivot_longer(c(logI, I, R, growth))

ggplot(df, aes(x = t, y = value)) +
  geom_line() +
  facet_wrap( ~ name, scales = "free_y", ncol = 1)

^{Created on 2024-07-02 with reprex v2.1.0}

[avallecam]

@avallecam can you check what it looks like when using rt = rt_opts(future = "project")?

Here it is

[seabbs]

The more concerning issue is that in the original question, the plot on the right has an axis up to 0.4 whilst on the left its 0.2 but looks visually the same in sample.

[avallecam]

The more concerning issue is that in the original question, the plot on the right has an axis up to 0.4 whilst on the left its 0.2 but looks visually the same in sample.

After updating to the new dist-interface, we made one misspecification as shown in this checklist box for the reporting_delay_fixed object. Should we expect this misspecification to get that effect? I found the previous versions of the same output from the original and current md-outputs.