epiforecasts / EpiNow2

Estimate Realtime Case Counts and Time-varying Epidemiological Parameters
https://epiforecasts.io/EpiNow2/dev/
Other
105 stars 31 forks source link

out of range doubling/halving central estimate from confidence intervals reported with default but not "vb" method #703

Open avallecam opened 2 weeks ago

avallecam commented 2 weeks ago

Summary: from the epinow() output, doubling/halving estimate, the range of confidence intervals does not include the central estimate.

Description: I would expect the central estimate to be included in the range. This is visible in the reference documentation of estimate_infection()

I think this could be method-specific, given that the issue appears with the default method, but not with the "vb" method.

Reproducible Steps:

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

# get example case counts
reported_cases <- example_confirmed[1:60]

# set an example generation time. In practice this should use an estimate
# from the literature or be estimated from data
generation_time <- Gamma(
  shape = Normal(1.3, 0.3),
  rate = Normal(0.37, 0.09),
  max = 14
)
# set an example incubation period. In practice this should use an estimate
# from the literature or be estimated from data
incubation_period <- LogNormal(
  meanlog = Normal(1.6, 0.06),
  sdlog = Normal(0.4, 0.07),
  max = 14
)
# set an example reporting delay. In practice this should use an estimate
# from the literature or be estimated from data
reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10)

# for more examples, see the "estimate_infections examples" vignette
def <- estimate_infections(reported_cases,
                           generation_time = generation_time_opts(generation_time),
                           delays = delay_opts(incubation_period + reporting_delay),
                           rt = rt_opts(prior = list(mean = 2, sd = 0.1)),
                           stan = stan_opts(samples = 1000, chains = 2)
)
#> Warning: Pareto k diagnostic value is 1.83. Resampling is disabled. Decreasing
#> tol_rel_obj may help if variational algorithm has terminated prematurely.
#> Otherwise consider using sampling instead.

# real time estimates
summary(def)
#>                             measure               estimate
#>                              <char>                 <char>
#> 1:           New infections per day     2309 (950 -- 5129)
#> 2: Expected change in daily reports      Likely decreasing
#> 3:       Effective reproduction no.      0.9 (0.59 -- 1.3)
#> 4:                   Rate of growth -0.03 (-0.15 -- 0.089)
#> 5:     Doubling/halving time (days)      -23 (7.8 -- -4.6)

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

EpiNow2 Version:

> packageVersion("EpiNow2")
[1] ‘1.5.2’

seabbs commented 1 week ago

I think this stems from the difficulty in defining the doubling time as it is a transformed quantity.

I believe it changes sign after passing through infinity and not 0 so -23 is in fact contained in these intervals.

I think this is probably a doc issue?

jamesmbaazam commented 6 days ago

I think this stems from the difficulty in defining the doubling time as it is a transformed quantity.

I believe it changes sign after passing through infinity and not 0 so -23 is in fact contained in these intervals.

I think this is probably a doc issue?