Initialisation error with data that has low initial incidence and little evidence of ongoing transmission.

[davidvilanova]

Hi,
Following the github exemple with my dataset i go this error,

I don´t know what is wrong ??

Below is my dataset that contains missing dates. Could the be the problem, should i include all days even if there is no data for that date ?

> reported_cases
           date confirm
  1: 2020-01-03       1
  2: 2020-01-21       1
  3: 2020-02-06       1
  4: 2020-02-24       1
  5: 2020-02-26       1
 ---                   
909: 2022-10-13      10
910: 2022-10-15       2
911: 2022-10-16      11
912: 2022-10-17      21
913: 2022-10-18       1

> 
estimates <- epinow(
  reported_cases = reported_cases,
  generation_time = generation_time,
  delays = delay_opts(incubation_period, reporting_delay),
  rt = rt_opts(prior = list(mean = 2, sd = 0.2)),
  stan = stan_opts(cores = 4),
  verbose = interactive(),
  zero_threshold  = 2
)

Logging threshold set at INFO for the EpiNow2 logger
Writing EpiNow2 logs to the console and: /tmp/RtmpVRQPCe/regional-epinow/2022-10-18.log
Logging threshold set at INFO for the EpiNow2.epinow logger
Writing EpiNow2.epinow logs to the console and: /tmp/RtmpVRQPCe/epinow/2022-10-18.log
DEBUG [2022-10-21 11:56:02] epinow: Running in exact mode for 2000 samples (across 4 chains each with a warm up of 250 iterations each) and 1034 time steps of which 7 are a forecast
starting worker pid=1338931 on localhost:11432 at 11:56:02.632
starting worker pid=1338945 on localhost:11432 at 11:56:02.922
starting worker pid=1338959 on localhost:11432 at 11:56:03.241
starting worker pid=1338978 on localhost:11432 at 11:56:03.532

SAMPLING FOR MODEL 'estimate_infections' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: normal_lpdf: Location parameter is -inf, but must be finite!  (in 'functions/rt.stan' at line 53; included from 'model_estimate_infections' at line 5)
  (in 'model_estimate_infections' at line 119)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: normal_lpdf: Location parameter is -inf, but must be finite!  (in 'functions/rt.stan' at line 53; included from 'model_estimate_infections' at line 5)
  (in 'model_estimate_infections' at line 119)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: normal_lpdf: Location parameter is -inf, but must be finite!  (in 'functions/rt.stan' at line 53; included from 'model_estimate_infections' at line 5)
  (in 'model_estimate_infections' at line 119)

.....
here are whatever error messages were returned
[[1]]
Stan model 'estimate_infections' does not contain samples.

[davidvilanova]

I have tried to add all missing days with a value of "0" but i still get the same error. Howver if add "1" to all newly added missing days it works. So my guess is that it fails when there is missing days or days with 0 infections

[seabbs]

Hi @davidvilanova,

Thanks for opening the issue. To make it easier to help could you provide a reprex (i.e a minimal example that reproduces the issue).

The model should be robust to the issues you describe. Either you have found a bug or there is an aspect of your data that doesn't match up to our model assumptions.

[davidvilanova]

Hi,
Sure here is a subset that fails for me.

reported_cases <- structure(list(date = structure(c(18264, 18282, 18298, 18316, 
18318, 18322, 18324, 18326, 18329, 18330, 18331, 18332, 18333, 
18334, 18335, 18336, 18337, 18338, 18339, 18340, 18341, 18342, 
18343, 18344, 18345, 18346, 18347, 18348, 18349, 18350, 18351, 
18352, 18353, 18354, 18355, 18356, 18357, 18358, 18359, 18360, 
18361, 18362, 18363, 18364, 18365, 18366, 18367, 18368, 18369, 
18370, 18371, 18372, 18373, 18374, 18375, 18376, 18377, 18378, 
18379, 18380), class = "Date"), confirm = c(1, 1, 1, 1, 1, 1, 
4, 1, 2, 1, 1, 2, 14, 12, 11, 28, 18, 38, 67, 64, 52, 20, 50, 
65, 63, 44, 49, 37, 63, 48, 67, 32, 25, 26, 59, 57, 56, 40, 56, 
20, 41, 38, 30, 17, 40, 44, 36, 31, 20, 12, 9, 15, 11, 14, 27, 
16, 3, 24, 23, 34)), row.names = c(NA, -60L), class = c("data.table", "data.frame"))

[seabbs]

Great thanks, and to be clear you are using the GitHub version of the package?

[davidvilanova]

I´m using EpiNow2_1.3.3,

thanks

[seabbs]

Great. If you have a chance posting a complete reprex (i.e data -> code -> error (I assume I can patch together code across posts here if not?)) along with session information would speed things up.

Do you also get this behaviour with 1.3.2?

[davidvilanova]

Ok here is the full code:

dataset:

reported_cases <- structure(list(date = structure(c(18264, 18282, 18298, 18316, 
18318, 18322, 18324, 18326, 18329, 18330, 18331, 18332, 18333, 
18334, 18335, 18336, 18337, 18338, 18339, 18340, 18341, 18342, 
18343, 18344, 18345, 18346, 18347, 18348, 18349, 18350, 18351, 
18352, 18353, 18354, 18355, 18356, 18357, 18358, 18359, 18360, 
18361, 18362, 18363, 18364, 18365, 18366, 18367, 18368, 18369, 
18370, 18371, 18372, 18373, 18374, 18375, 18376, 18377, 18378, 
18379, 18380), class = "Date"), confirm = c(1, 1, 1, 1, 1, 1, 
4, 1, 2, 1, 1, 2, 14, 12, 11, 28, 18, 38, 67, 64, 52, 20, 50, 
65, 63, 44, 49, 37, 63, 48, 67, 32, 25, 26, 59, 57, 56, 40, 56, 
20, 41, 38, 30, 17, 40, 44, 36, 31, 20, 12, 9, 15, 11, 14, 27, 
16, 3, 24, 23, 34)), row.names = c(NA, -60L), class = c("data.table", "data.frame"))

Code:

generation_time <- get_generation_time(
  disease = "SARS-CoV-2", source = "ganyani"
)
incubation_period <- get_incubation_period(
  disease = "SARS-CoV-2", source = "lauer"
)
reporting_delay <- list(
  mean = convert_to_logmean(2, 1), mean_sd = 0.1,
  sd = convert_to_logsd(2, 1), sd_sd = 0.1,
  max = 10
)

estimates <- epinow(
  reported_cases = reported_cases,
  generation_time = generation_time,
  delays = delay_opts(incubation_period, reporting_delay),
  rt = rt_opts(prior = list(mean = 2, sd = 0.2)),
  stan = stan_opts(cores = 4),
  verbose = interactive()
  )

Partial error code

Chain 3: Rejecting initial value:
Chain 3:   Error evaluating the log probability at the initial value.
Chain 3: Exception: Exception: normal_lpdf: Location parameter is -inf, but must be finite!  (in 'functions/rt.stan' at line 53; included from 'model_estimate_infections' at line 5)
  (in 'model_estimate_infections' at line 119)

Chain 3: Rejecting initial value:
Chain 3:   Error evaluating the log probability at the initial value.
Chain 3: Exception: Exception: normal_lpdf: Location parameter is -inf, but must be finite!  (in 'functions/rt.stan' at line 53; included from 'model_estimate_infections' at line 5)
  (in 'model_estimate_infections' at line 119)

[davidvilanova]

i didn´t try with 1.3.2

[davidvilanova]

> sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.8.0 EpiNow2_1.3.3   stringr_1.4.1   httr_1.4.4      elastic_1.2.0   ggplot2_3.3.6  
 [7] dplyr_1.0.9     EpiEstim_2.2-4  tidyr_1.2.0     epitrix_0.2.2   incidence_1.7.3 outbreaks_1.9.0
oaded via a namespace (and not attached):
 [1] jsonlite_1.8.0        splines_4.2.1         R.utils_2.12.0        StanHeaders_2.21.0-7 
 [5] RcppParallel_5.1.5    assertthat_0.2.1      triebeard_0.3.0       urltools_1.7.3       
 [9] stats4_4.2.1          globals_0.16.1        pillar_1.8.1          lattice_0.20-45      
[13] quantreg_5.94         glue_1.6.2            digest_0.6.29         runner_0.4.2         
[17] colorspace_2.0-3      Matrix_1.5-1          R.oo_1.25.0           plyr_1.8.7           
[21] pkgconfig_2.0.3       httpcode_0.3.0        rstan_2.21.7          SparseM_1.81         
[25] listenv_0.8.0         purrr_0.3.4           coarseDataTools_0.6-6 patchwork_1.1.2      
[29] scales_1.2.1          processx_3.7.0        MatrixModels_0.5-0    tibble_3.1.8         
[33] generics_0.1.3        ellipsis_0.3.2        withr_2.5.0           cli_3.3.0            
[37] survival_3.4-0        magrittr_2.0.3        crayon_1.5.1          mcmc_0.9-7           
[41] ps_1.7.1              R.methodsS3_1.8.2     future_1.28.0         fansi_1.0.3          
[45] parallelly_1.32.1     MASS_7.3-58           truncnorm_1.0-8       pkgbuild_1.3.1       
[49] progressr_0.11.0      loo_2.5.1             prettyunits_1.1.1     tools_4.2.1          
[53] fitdistrplus_1.1-8    data.table_1.14.2     matrixStats_0.62.0    formatR_1.12         
[57] lifecycle_1.0.1       MCMCpack_1.6-3        munsell_0.5.0         lambda.r_1.2.4       
[61] callr_3.7.2           compiler_4.2.1        rlang_1.0.4           futile.logger_1.4.3  
[65] grid_4.2.1            rstudioapi_0.14       gtable_0.3.0          codetools_0.2-18     
[69] inline_0.3.19         DBI_1.1.3             curl_4.3.2            reshape2_1.4.4       
[73] R6_2.5.1              gridExtra_2.3         rstantools_2.2.0      future.apply_1.9.1   
[77] utf8_1.2.2            futile.options_1.0.1  stringi_1.7.8         parallel_4.2.1       
[81] crul_1.2.0            Rcpp_1.0.9            vctrs_0.4.1           tidyselect_1.1.2     
[85] coda_0.19-4          
> 

[seabbs]

So the core issue is that the first part of your data starts with 3 months with effectively zero cases (i.e due to imports and not local transmission). As you will see from our documentation, we are using a renewal equation, so previous cases generate new cases in our model. Here there is very little or no transmission, so the model is unstable.

I think the error comes from the initialisation of the initial conditions. Ideally this would be handled more robustly for these settings. Even if this is handled however I think estimates from the current model would be very unstable just due to the generating process of the data being far from our assumed process.

If you filter to the start of more transmission-like behaviour (i.e the 1st of March) I see successful model fits.

Obviously, your setting is one we might see at the start of an outbreak. To support it we would need a constant importation factor (or treat this as known data). Both of these are options (with the latter being what is done in EpiEstim and the former hhh4). I've explored fitting an importation factor previously but unfortunately in most setting, it is not identifiable.

If the filtering approach is not sufficient and you'd like guidance on alternative solutions (and potentially adding these to the package) happy to provide guidance.

[davidvilanova]

Ok, great, i see !!!
It works with date > 1-March. I still however find it extremely slow to run with a dataframe of 912 rows.

[seabbs]

Yes - that is because this is a method that is keyed for real-time usage and not retrospective usage by default. If you want retrospective estimates, and it sounds like you do if trying to fit to that much data, I'd explore some of the faster options (like changing the mean function of the gaussian process in rt_opts().

[davidvilanova]

Great Thanks !!
just notice that of you adjuts the cores number it is not taking this into account?. For example if i switch to 20 cores i get. Should be 20 chains?

That´s the output with 20 cores:
DEBUG [2022-10-24 13:34:59] epinow: Running in exact mode for 1000 samples **(across 4 chains** each with a warm up of 250 iterations each) and 980 time steps of which 7 are a forecast

[seabbs]

you'd need to also change the number of MCMC chains to match. It isn't ideal as all chains will need to warm-up (i.e so won't give a linear speed up) but will definitely help get a richer sample of the posterior in less time.