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

[avallecam]

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]

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]

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?

[sbfnk]

The doubling time is defined in

EpiNow2/R/report.R

Line 60 in fc92e60

doubling_time <- function(r) {

It's calculated from the growth rate $r$ as $\frac{\log 2}{r}$. If $r$ is estimated as (-0.15, 0.089) then for the doubling time you go from $\log(2)/(-0.15) = -4.5$ to $-\infty$ as $r$ approaches 0 (via the central estimate $\log(2)/-0.03=-23.1$) and then as $r$ crosses 0 you come back from $+\infty$ down to $\log(2)/0.089 = 7.8.

So what we're saying is that the doubling time is either a halving time (i.e. it's negative) of more than 4.5 days or a doubling time of more than 7.8 days, with a central estimate being a halving time of 23 days (if growth is 0 then the doubling/halving time is infinite).

library("ggplot2")

r <- seq(-0.15, 0.089, by = 0.01)
dt <- log(2) / r

df <- data.frame(r = r, dt = dt)

ggplot(df, aes(x = r, y = dt)) +
  geom_line() +
  geom_vline(xintercept = -0.03, linetype = "dashed")

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