Skip to content

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
avallecam opened this issue Jul 2, 2024 · 6 comments

Comments

@avallecam
Copy link

avallecam commented Jul 2, 2024

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

Description:

image (4)

  • 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
Copy link
Contributor

seabbs commented Jul 2, 2024

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

@sbfnk
Copy link
Contributor

sbfnk commented Jul 2, 2024

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
Copy link
Contributor

sbfnk commented Jul 2, 2024

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)

image

Created on 2024-07-02 with reprex v2.1.0

@avallecam
Copy link
Author

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

Here it is
image

@seabbs
Copy link
Contributor

seabbs commented Jul 3, 2024

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
Copy link
Author

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.

@epiforecasts epiforecasts locked and limited conversation to collaborators Jul 4, 2024
@seabbs seabbs converted this issue into discussion #711 Jul 4, 2024

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants