Hi Anne Cori,

I am using EpiEstim to estimate the instantaneous (case) reproduction number for 2009 pandemic influenza A (H1N1) in mainland China. The following are my code:

rm(list = ls())

load(url("http://tonytsai.name/confirmed_pdm_dec.rda"))

# instantaneous reproduction number estimation for pandemic --------------------
# using ParametricSI method
# the instantaneous reproduction number can be estimated after May 22nd, 2009
EstimateR(dec$cases, T.Start = 22:359, T.End = 28:365, method = "ParametricSI", 
          Mean.SI = 2.6, Std.SI = 1.3, plot = TRUE, leg.pos = xy.coords(1, 3))
# case reproduction number estimaion for pandemic ------------------------------
# using the Wallinga and Teunis method
WT(dec$cases, T.Start = 20:100, T.End = 26:106, method = "ParametricSI", Mean.SI = 2.6, 
   Std.SI = 1.3, plot = TRUE, nSim = 100)

The instantaneous reproduction number can be successfully estimated, but the WT function failed and the R session aborted.

Could you help me to fix the problem with WT function? Thank you very much.

There are two branches that are ahead of master, new-version and release. It is confusing why both of these should be ahead of master. When comparing these, it appears that release may be slightly ahead of new-version and should be favored: https://github.com/annecori/EpiEstim/compare/new-version..release

proposed changes to allow a Date column to be specified in I, which is then used for plotting --> addresses issue #12

also, added errors when the estimation is performed to early or too late --> addresses issue #15 and #19

finally, also allowed EstimateR and WT to take incidence objects (from class incidence from package incidence) as arguments --> addresses issue #13

Dr. Robin Thomas asked me to submit this bug report. There is an error in the EpiEstim app which causes the 75% & 25% confidence intervals to show as identical.

t_start | t_end | Mean(R) | Std(R) | Quantile.0.025(R) | Quantile.0.05(R) | Quantile.0.25(R) | Median(R) | Quantile.0.75(R) | Quantile.0.95(R) | Quantile.0.975(R) -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- 2 | 8 | 1.676272 | 0.976486 | 0.338931 | 0.449985 | 3.524117 | 1.492907 | 3.524117 | 3.524117 | 4.034139 3 | 9 | 2.584599 | 1.163681 | 0.839038 | 1.020117 | 4.750958 | 2.408954 | 4.750958 | 4.750958 | 5.33603 4 | 10 | 2.940135 | 1.134309 | 1.165467 | 1.355705 | 5.003363 | 2.791074 | 5.003363 | 5.003363 | 5.543205 5 | 11 | 2.29199 | 0.892125 | 0.902246 | 1.056114 | 3.929101 | 2.176504 | 3.929101 | 3.929101 | 4.352027 6 | 12 | 2.222825 | 0.803635 | 0.943369 | 1.096749 | 3.693079 | 2.123938 | 3.693079 | 3.693079 | 4.045335 7 | 13 | 2.13141 | 0.730564 | 0.960056 | 1.099175 | 3.455796 | 2.046818 | 3.455796 | 3.455796 | 3.792869 8 | 14 | 3.563511 | 0.88762 | 2.06355 | 2.251539 | 5.157335 | 3.483487 | 5.157335 | 5.157335 | 5.525408 9 | 15 | 2.845731 | 0.687544 | 1.678171 | 1.830796 | 4.07437 | 2.7868 | 4.07437 | 4.07437 | 4.353653 10 | 16 | 2.918638 | 0.626415 | 1.850218 | 1.98553 | 4.028234 | 2.865315 | 4.028234 | 4.028234 | 4.293019

Lauren McGough (@unrealmcg) and I have been doing some simple tests to compare Rt methods on synthetic data. We've been running into errors when calling the wallinga_teunis() function in EpiEstim, of the form values must be length <A,> but FUN(X[[1]]) result is length <B>.

This only happens when n_sim > 0. If n_sim == 0—skipping the CIs—it seems to be fine.

E.g.:

Error in vapply(seq_len(config$n_sim), function(i) draw_one_set_of_ancestries(),  : 
  values must be length 19889,
 but FUN(X[[1]]) result is length 19885
Calls: wallinga_teunis -> t -> vapply
Execution halted

That error came from the following code, with inline data (just generated from a stochastic SEIR model):

library(EpiEstim)

incidence <- c(
  1, 3, 2, 2, 2, 1, 1, 1, 1, 1, 4, 1, 4, 3, 2, 2, 2, 3, 7, 8, 3, 0, 1, 0, 3, 3, 3, 2, 1, 1, 3, 1, 3, 2, 0, 0, 3, 2, 0, 1, 2, 0, 2, 2, 1, 1, 2, 1, 2, 2, 1, 1, 2, 3, 5, 5, 5, 3, 4, 5, 3, 6, 2, 3, 10, 8, 7, 7, 11, 5, 7, 11, 7, 4, 12, 10, 9, 13, 10, 12, 9, 5, 8, 9, 6, 8, 11, 9, 12, 12, 7, 12, 9, 15, 10, 8, 13, 13, 19, 8, 5, 14, 15, 10, 15, 12, 17, 14, 13, 13, 14, 16, 16, 14, 11, 13, 19, 21, 15, 15, 20, 14, 11, 23, 12, 20, 21, 18, 18, 19, 18, 20, 20, 17, 18, 31, 28, 13, 29, 20, 24, 31, 25, 29, 23, 33, 24, 27, 30, 26, 26, 24, 25, 21, 28, 41, 31, 32, 47, 29, 37, 36, 35, 35, 35, 46, 41, 37, 38, 28, 41, 35, 35, 38, 20, 31, 38, 42, 35, 31, 42, 39, 47, 30, 57, 33, 40, 29, 28, 41, 34, 33, 42, 48, 32, 38, 33, 46, 45, 41, 42, 46, 42, 39, 52, 43, 46, 44, 33, 45, 56, 36, 54, 51, 52, 45, 51, 57, 55, 59, 60, 45, 46, 56, 37, 49, 58, 38, 55, 47, 60, 51, 41, 51, 36, 63, 35, 43, 57, 60, 43, 60, 60, 51, 44, 51, 64, 65, 75, 68, 65, 66, 62, 69, 57, 67, 67, 69, 68, 78, 60, 72, 64, 66, 61, 67, 55, 71, 82, 60, 78, 77, 70, 76, 64, 63, 62, 58, 75, 80, 71, 80, 67, 57, 67, 63, 81, 77, 77, 72, 74, 69, 64, 83, 66, 77, 73, 62, 64, 82, 72, 72, 58, 56, 66, 86, 68, 70, 63, 71, 60, 61, 57, 54, 54, 60, 58, 60, 62, 68, 46, 70, 75, 59, 73, 58, 67, 50, 66, 59, 69, 68, 63, 76, 62, 62, 58, 66, 60, 75, 60, 78, 63, 53, 70, 66, 71, 46, 61, 66, 72, 75, 83, 64, 73, 64, 55, 88, 63, 66, 67, 66, 78, 62, 71, 70, 77, 65, 45, 76, 73, 72, 53, 50, 68, 65, 66, 44, 52, 59, 77, 52, 66, 61, 66, 64, 68, 59, 64, 51, 46, 57, 61, 52, 44, 58, 48, 40, 48, 55, 62, 42, 50, 53, 39, 53, 50, 49, 53, 49, 43, 44, 49, 44, 43, 42, 39, 37, 37, 34, 41, 50, 46, 30, 43, 45, 35, 27, 37, 45, 32, 46, 26, 26, 32, 27, 34, 34, 23, 33, 36, 28, 36, 33, 32, 29, 38, 31, 30, 30, 38, 27, 34, 38, 34, 19, 27, 35, 32, 28, 36, 26, 25, 33, 23, 26, 28, 20, 27, 24, 25, 20, 28, 21, 20, 26, 24, 19, 16, 21, 22, 17, 23, 22, 17, 24, 30, 17, 16, 18, 16, 15, 17, 18, 16, 14, 18, 21, 18, 14, 19, 17, 17, 10, 19, 19, 14, 13, 15, 9, 9, 10, 13, 10, 9, 13, 8, 10, 14, 9, 9, 10, 5, 17, 14, 10, 14, 14, 5, 15, 12, 9, 11, 18, 12, 11, 12, 14, 13, 13, 10, 10, 17, 15, 7, 13, 11, 8, 7, 9, 9, 7, 9, 6, 10, 14, 10, 7, 3, 5, 11, 9, 4, 7, 5, 5, 7, 5, 9, 8, 6, 3, 4, 8, 6, 6, 8, 5, 5, 5, 6, 8, 4, 3, 7, 8, 7, 3, 5, 7, 7, 4, 2, 4, 7, 1, 2, 3, 3, 5, 4, 3, 2, 4, 5, 1, 3, 1, 3, 1, 3, 3, 4, 2, 6, 0, 2, 6, 7, 4, 4, 4, 2, 0, 6, 0, 1, 2, 3, 0, 1, 2, 5, 3, 5, 3, 1, 1, 3, 1, 3, 1, 4, 2, 4, 3, 2, 2, 3, 3, 1, 1, 3, 6, 3, 2, 1, 2, 3, 4, 3, 2, 0, 2, 4, 3, 4, 0, 5, 2, 1, 1, 4, 1, 1, 2, 2, 5, 2, 1, 1, 4, 1, 3, 3, 4, 3, 5, 3, 3, 5, 4, 2, 0, 2, 3, 5, 3, 2, 7, 1, 1, 2, 1, 2, 1, 1, 3, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 2, 0, 0, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 0, 0, 2, 1, 1, 1, 0, 0, 2, 2, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 3, 1, 3, 1, 1, 1, 2, 0, 3, 1, 0, 2, 1, 1, 0, 0, 2, 0, 1, 1, 0, 0, 2, 1, 0, 2, 1, 0, 0, 0, 0, 2, 1, 1, 0, 0, 1, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0
)
n_t <- length(incidence)

mean_serial_int <- 1/(1.1 / 5) + 3
std_serial_int <- mean_serial_int

window_size <- 11

t_start <- seq(2, n_t - 20)
t_end <- t_start + window_size - 1
wt_result <- wallinga_teunis(
  incidence, method="parametric_si",
  config = list(
    t_start = t_start,
    t_end = t_end,
    mean_si = mean_serial_int,
    std_si = std_serial_int,
    n_sim = 10
  )
)

library(EpiEstim)
data(Flu2009)
T <- nrow(Flu2009$incidence)
t_start <- seq(2, T-6) # starting at 2 as conditional on the past observations
t_end <- t_start + 6 # adding 6 to get 7-day windows as bounds included in window
res_weekly <- EpiEstim::estimate_R(Flu2009$incidence, 
                         method="parametric_si",
                         config = EpiEstim::make_config(list(
                             t_start = t_start,
                             t_end = t_end,
                             mean_si = 2.6, 
                             std_si = 1.5)))
res_weekly$R

results in :

t_start t_end   Mean(R)     Std(R) Quantile.0.025(R) Quantile.0.05(R) Quantile.0.25(R) Median(R) Quantile.0.75(R)
1        2     8 1.7357977 0.40913143        1.02874370       1.12193325        2.4589724 1.7037612        2.4589724
2        3     9 1.7491678 0.36472669        1.10882231       1.19547993        2.3891206 1.7238839        2.3891206

Other quantiles look OK

Related to #40,

The new version of EpiEstim currently has no tests and that's.... not good. In fact, with the current master branch, Example 2 fails.

Regarding tests, the current setup is relatively reasonable since they do not rely on randomization to generate the data, but we need to find out why Example 2 is no good.

This could have been caught earlier with continuous integration, so I would suggest to use the following to create it.

usethis::use_travis()
usethis::use_appveyor()

Hi Anne Cori,

I am using EpiEstim to estimate the instantaneous (case) reproduction number during post-pandemic period for 2009 pandemic influenza A (H1N1) in mainland China. The EstimateR function successfully estimated the R(t); however the maximal estimation of R(t) is 47.5, which is so large that I don't think it makes sense. Could you help me to explain why such a large estimation of R(t) could be produced? Thank you very much.

> rm(list = ls())
> 
> load(url("http://tonytsai.name/confirmed_post-pdm_dec.rda"))
> 
> # instantaneous reproduction number estimation for post-pandemic --------------------
> # using ParametricSI method
> # the instantaneous reproduction number can be estimated after May 2nd, 2010
> x <- EstimateR(dec$cases, T.Start = 2:359, T.End = 8:365, method = "ParametricSI", 
+                Mean.SI = 2.6, Std.SI = 1.3, plot = TRUE, leg.pos = xy.coords(1, 3))
> max(x$R$`Mean(R)`)
[1] 47.54329

Pull Request Closes #92

• Fixes a bug where draw_one_set_of_ancestries would return a result of the wrong length. It would calculate the length based on the time window, but everything else is based on T. I am not familiar with the actual maths involved here, so please do check this is correct.

• Fixes a bug where ot was not defined.

How has this been tested Examples were given in #92, and these now work correctly.

Checklist

• [X] I have added tests to prove my changes work
• [X] I have added documentation where required
• [X] I have updated NEWS.md with a short description of my change

This will fix #70, but I've modified it so that it takes from a specific time window of R:

library("EpiEstim")
#> Registered S3 methods overwritten by 'ggplot2':
#>   method         from 
#>   [.quosures     rlang
#>   c.quosures     rlang
#>   print.quosures rlang
data("Flu2009")

res <- estimate_R(incid = Flu2009$incidence, 
                  method = "non_parametric_si",
                  config = make_config(list(si_distr = Flu2009$si_distr)))
#> Default config will estimate R on weekly sliding windows.
#>     To change this change the t_start and t_end arguments.

hist(sample_posterior_R(res, n = 5000, window = 1L), col = "grey",
     main = "5000 samples of R from the first weekly window",
     xlab = "R",
     xlim = c(0, 4))

hist(sample_posterior_R(res, n = 5000, window = 10L), col = "grey",
     main = "5000 samples of R from the tenth weekly window",
     xlab = "R",
     xlim = c(0, 4))

win_col <- ifelse(seq(nrow(res$R)) %in% c(1, 10), "red", "black")
plot(res, "R") + ggplot2::geom_point(color = win_col)

^{Created on 2019-06-06 by the reprex package (v0.3.0)}

Sorry I've been absent on this. It would be good to tag the new version as it was released to CRAN. I think tagging the most recent commit with 2.2-3 would be sufficient.

In line with https://github.com/mrc-ide/EpiEstim/issues/152, it would be useful to provide an S3 method for incidence2 inputs. The incidence2 package is meant as a replacement for incidence, and offers more flexibility. Some issues to think about / handle:

• handle multiple stratifications
• handle non-days time intervals (may need postponing into a separate issue)

Turning the main function into a generic will facilitate providing dedicated functions for different types of inputs, e.g. an integer vector, and incidence, or an incidence2 object.

Hi!

First of all thank you so much for this great package! I downloaded EpiEstim a few months ago using install.packages() and I've only been using the estimate_R function so far and that has worked fine. Today, I needed to use the estimate_advantage function but that gave me an error saying that the function couldn't be found. I couldn't access the vignette associated with it (MV_EpiEstim_vignette) either. I tried uninstalling and reinstalling it but that didn't fix the problem so I uninstalled it again and then installed it using devtools::install_github instead and that worked. I'm not sure if I did something weird when I installed it initially, but I thought I should let you know!

Best, Anjalika

estimate_R currently returns the mean and standard deviation of R, which then can be used to draw samples from the Gamma. It would be convenient to have an option that the posterior draws from estimate_R are returned directly for subsequent use in the projections package